1module lamSpectralTransform_mod
2 ! MODULE lamSpectralTransform_mod (prefix='lst' category='4. Data Object transformations')
3 !
4 !:Purpose: Bi-Fourier spectral transform for limited-area applications.
5 ! Depends on ffft8 and setfft8 routines in ARMNLIB.
6 !
7 use mpi
8 use midasMpi_mod
9 use MathPhysConstants_mod
10 use earthConstants_mod
11 use utilities_mod
12 implicit none
13 save
14 private
15
16 ! public derived type
17 public :: struct_lst
18 ! public procedures
19 public :: lst_Setup, lst_Laplacian, lst_VarTransform
20 public :: lst_ReshapeTrunc ! only for standalone tests
21
22 type :: struct_lst
23 real(8), allocatable :: k_r8(:) ! (real) Total Wavenumber associated with each
24 ! nla spectral coefficient
25 integer, allocatable :: k(:) ! Total Wavenumber bin associated with each
26 ! nla spectral coefficient
27 integer, allocatable :: m(:) ! Wavenumber in x associated with each
28 ! nla spectral coefficient
29 integer, allocatable :: n(:) ! Wavenumber in y associated with each
30 ! nla spectral coefficient
31 real(8), allocatable :: Weight(:) ! Weight associated with each
32 ! nla spectral coefficient
33 integer, allocatable :: nePerK(:) ! Number of spectral element in each
34 ! total wavenumber bands
35 integer, allocatable :: nePerKglobal(:) ! Number of spectral element in each
36 ! total wavenumber bands over ALL PROCESSORS
37 integer, allocatable :: ilaFromEK(:,:) ! ila index associated to each spectral element
38 ! of total wavenumber band
39 real(8), allocatable :: NormFactor(:,:)
40 real(8), allocatable :: NormFactorAd(:,:)
41 integer :: nlaGlobal ! First dimension of VAR global spectral array
42 integer, allocatable :: ilaGlobal(:) ! Position of the local vector element in the global
43 ! vector
44 integer :: nphase ! Second dimension of VAR spectral array
45 integer :: ni
46 integer :: nj
47 integer :: nip
48 integer :: njp
49 integer :: mmax, nmax
50 integer :: ktrunc
51 integer :: nla ! First dimension of VAR spectral array
52 integer :: maxnla
53 integer :: mymBeg, mymEnd, mymSkip, mymCount, mymActiveCount,maxmActiveCount
54 integer :: mynBeg, mynEnd, mynSkip, mynCount
55 integer, allocatable :: nla_Index(:,:)
56 real(8), allocatable :: lapxy(:)
57 real(8), allocatable :: ilapxy(:)
58 integer, allocatable :: allmBeg(:), allmEnd(:), allmSkip(:)
59 integer, allocatable :: mymIndex(:)
60 integer, allocatable :: allnBeg(:), allnEnd(:), allnSkip(:)
61 integer, allocatable :: mynIndex(:)
62 logical :: allocated = .false.
63 integer :: latPerPE, latPerPEmax, myLatBeg, myLatEnd
64 integer :: lonPerPE, lonPerPEmax, myLonBeg, myLonEnd
65 integer :: myLevBeg, myLevEnd, myLevCount, maxLevCount
66 integer, allocatable :: allLatBeg(:), allLatEnd(:), allLatPerPE(:)
67 integer, allocatable :: allLonBeg(:), allLonEnd(:), allLonPerPE(:)
68 integer, allocatable :: allLevBeg(:), allLevEnd(:)
69 integer, allocatable :: KfromMNglb(:,:)
70 character(len=10) :: MpiMode
71 character(len=3) :: gridDataOrder ! Ordering the gridded data: 'ijk' or 'kij'
72 integer :: sendType_LevToLon, recvType_LevToLon
73 integer :: sendType_LonToLev, recvType_LonToLev
74 logical :: lonLatDivisible
75 end type struct_lst
76
77 ! TransformType = 'SinCos'
78 integer, parameter :: nip = 2 ! Padding
79 integer, parameter :: njp = 2 ! Padding
80 integer, parameter :: nphase = 4 ! For Sin&Cos we have a) sin(m)sin(n)
81 ! b) cos(m)cos(n)
82 ! c) sin(m)cos(n)
83 ! d) cos(m)sin(n)
84
85 logical, parameter :: verbose = .false.
86
87contains
88
89 !--------------------------------------------------------------------------
90 ! lst_setup
91 !--------------------------------------------------------------------------
92 subroutine lst_Setup(lst, ni_in, nj_in, dlon_in, ktrunc_in, &
93 MpiMode, maxlevels_opt, gridDataOrder_opt)
94 implicit none
95
96 ! Arguments:
97 type(struct_lst), intent(out) :: lst ! Parameters available to the outside world
98 integer, intent(in) :: ni_in ! Global Grid point data horizontal dimensions
99 integer, intent(in) :: nj_in ! Global Grid point data horizontal dimensions
100 character(len=*), intent(in) :: MpiMode ! MPI Strategy
101 real(8), intent(in) :: dlon_in ! Grid Spacing in Radians
102 integer, intent(in) :: ktrunc_in ! Spectral Truncation (global)
103 integer, optional, intent(in) :: maxlevels_opt ! Number of levels; Only needed when MpiMode = LatLev
104 character(len=*), optional, intent(in) :: gridDataOrder_opt ! 'ijk' or 'kij'
105
106 ! Locals:
107 real(8), allocatable :: Kr8fromMN(:,:)
108 integer, allocatable :: KfromMN(:,:)
109 integer, allocatable :: my_KfromMNglb(:,:)
110 integer :: kref, mref, nref
111 integer :: m, n, k, kMax, ila, nfact_lon, nfact_lat
112 integer :: ier, ilaglb, i, j, p
113 real(8) :: dlon, dx2, fac, ca, cp, cb, cq, r
114 real(8) :: NormFactor1, NormFactor2, NormFactor3
115 real(8) :: NormFactorAd1, NormFactorAd2, NormFactorAd3
116 real(8) :: factor, factorAd
117 character(len=60) :: kreftype
118 logical :: divisibleLon, divisibleLat
119 integer(kind=MPI_ADDRESS_KIND) :: lowerBound, extent
120 integer :: realSize, sendType, recvType, ierr
121
122 !
123 !- 1. Set variables needed by the LAM Spectral Transform in VAR
124 !
125 if (verbose) write(*,*) 'Entering lst_Setup'
126
127 if (lst%allocated) then
128 call utl_abort('lst_setup: this structure is already allocated')
129 end if
130
131 kreftype = 'MAX' ! hardwired
132
133 lst%ni = ni_in
134 lst%nj = nj_in
135 lst%nphase = nphase
136 lst%nip = nip
137 lst%njp = njp
138
139 ! 1.1 Check grid dimensions and set padding for the RPN DFT routines
140
141 ! We need to padd the input array such as ...
142 ! O O O O O O O O O
143 ! O O O O O O O O O
144 ! X X X X X X X O O
145 ! X X X X X X X O O
146 ! X X X X X X X O O
147 ! X X X X X X X O O
148
149 if (mod(lst%ni,2) /= 0 .or. mod(lst%nj,2) /= 0) then
150 write(*,*) ' The regular Sin & Cos Transform requires that', &
151 ' dimensions be EVEN. Fields MUST be periodic' , &
152 ' but the last colum and row MUST NOT BE a ' , &
153 ' repetition of the first colum and row. '
154 call utl_abort('lst_setup')
155 end if
156
157 nfact_lon = lst%ni
158 call ngfft(nfact_lon) ! INOUT
159 nfact_lat = lst%nj
160 call ngfft(nfact_lat) ! INOUT
161
162 if (nfact_lon /= lst%ni .or. nfact_lat /= lst%nj) then
163 if (nfact_lon /= lst%ni) then
164 write(*,*) 'Error: A fast transform cannot be used in X'
165 write(6,6130) lst%ni, nfact_lon
166 end if
167 if (nfact_lat /= lst%nj) then
168 write(*,*) 'Error: A fast transform cannot be used in Y'
169 write(6,6140) lst%nj, nfact_lat
170 end if
171 call utl_abort('lst_setup')
172 end if
173
1746130 FORMAT('N = ni = ', I4,' the nearest factorizable N = ',I4)
1756140 FORMAT('N = nj = ', I4,' the nearest factorizable N = ',I4)
176
177 !- 1.2 Maximum of integer wavenumbers in x and y directions
178 lst%mmax = lst%ni/2
179 lst%nmax = lst%nj/2
180
181 write(*,'(A,f8.1)') ' lst_Setup: Your grid spacing (in km) = ', ec_ra*dlon_in/1000.0
182 write(*,*) ' Max wavenumbers in x-axis = ', lst%mmax
183 write(*,*) ' Max wavenumbers in y-axis = ', lst%nmax
184
185 !- 1.3 MPI Strategy
186
187 lst%MpiMode = MpiMode
188 select case (trim(lst%MpiMode))
189 case ('NoMpi')
190 !- 1.3.1 No MPI
191
192 ! range of LONS handled by this processor (ALL) in GRIDPOINT SPACE
193 lst%lonPerPE = lst%ni
194 lst%myLonBeg = 1
195 lst%myLonEnd = lst%ni
196
197 ! range of LATS handled by this processor (ALL) in GRIDPOINT SPACE
198 lst%latPerPE = lst%nj
199 lst%myLatBeg = 1
200 lst%myLatEnd = lst%nj
201
202 ! range of M handled by this processor (ALL) in SPECTRAL SPACE
203 lst%mymBeg = 0
204 lst%mymEnd = lst%mmax
205 lst%mymSkip = 1
206 lst%mymCount = lst%mmax + 1
207
208 ! range of N handled by this processor (ALL) in SPECTRAL SPACE
209 lst%mynBeg = 0
210 lst%mynEnd = lst%nmax
211 lst%mynSkip = 1
212 lst%mynCount = lst%nmax + 1
213
214 ! set a dummy range of LEVELS handled by this processor
215 lst%myLevBeg = -1
216 lst%myLevEnd = -1
217 lst%myLevCount = 0
218
219 case ('LatLonMN')
220 !- 1.3.2 MPI 2D: Distribution of lon/lat tiles (gridpoint space) and n/m (spectral space)
221
222 ! range of LONS handled by this processor in GRIDPOINT SPACE
223 call mmpi_setup_lonbands(lst%ni, & ! IN
224 lst%lonPerPE, lst%lonPerPEmax, & ! OUT
225 lst%myLonBeg, lst%myLonEnd, & ! OUT
226 divisible_opt=divisibleLon) ! OUT
227
228 ! range of LATS handled by this processor in GRIDPOINT SPACE
229 call mmpi_setup_latbands(lst%nj, & ! IN
230 lst%latPerPE, lst%latPerPEmax, & ! OUT
231 lst%myLatBeg, lst%myLatEnd, & ! OUT
232 divisible_opt=divisibleLat) ! OUT
233
234 lst%lonLatDivisible = (divisibleLon .and. divisibleLat)
235 if(mmpi_myid == 0) write(*,*) 'lst_setup: lonLatDivisible = ', lst%lonLatDivisible
236
237 ! range of M handled by this processor in SPECTRAL SPACE
238 call mmpi_setup_m(lst%mmax, & ! IN
239 lst%mymBeg, lst%mymEnd, lst%mymSkip, lst%mymCount) ! OUT
240
241 ! range of N handled by this processor in SPECTRAL SPACE
242 call mmpi_setup_n(lst%nmax, & ! IN
243 lst%mynBeg, lst%mynEnd, lst%mynSkip, lst%mynCount) ! OUT
244
245 ! range of LEVELS TEMPORARILY handled by this processor DURING THE SPECTRAL TRANSFORM
246 if (.not.present(maxlevels_opt)) then
247 call utl_abort('lst_setup: ERROR, number of levels must be specified with MpiMode LatLonMN')
248 end if
249 ! 2D MPI decomposition: split levels across npex
250 call mmpi_setup_levels(maxlevels_opt, & ! IN
251 lst%myLevBeg,lst%myLevEnd,lst%myLevCount) ! OUT
252
253 case default
254 write(*,*)
255 write(*,*) 'Error: MpiMode Unknown ', trim(MpiMode)
256 call utl_abort('lst_setup')
257 end select
258
259 write(*,*)
260 write(*,*) ' I am processor ', mmpi_myid+1, ' on a total of ', mmpi_nprocs
261 write(*,*) ' mband info = ', lst%mymBeg, lst%mymEnd, lst%mymSkip, lst%mymCount
262 write(*,*) ' nband info = ', lst%mynBeg, lst%mynEnd, lst%mynSkip, lst%mynCount
263 write(*,*) ' level info = ', lst%myLevBeg, lst%myLevEnd, lst%myLevCount
264
265 ! Set M index
266 allocate(lst%mymIndex(lst%mymBeg:lst%mymEnd))
267 lst%mymIndex(:)=0
268 do m = lst%mymBeg, lst%mymEnd, lst%mymSkip
269 if (m == lst%mymBeg) then
270 lst%mymIndex(m) = 1
271 else
272 lst%mymIndex(m) = lst%mymIndex(m-lst%mymSkip) + 1
273 end if
274 end do
275
276 ! Set N index
277 allocate(lst%mynIndex(lst%mynBeg:lst%mynEnd))
278 lst%mynIndex(:)=0
279 do n = lst%mynBeg, lst%mynEnd, lst%mynSkip
280 if (n == lst%mynBeg) then
281 lst%mynIndex(n) = 1
282 else
283 lst%mynIndex(n) = lst%mynIndex(n-lst%mynSkip) + 1
284 end if
285 end do
286
287 if (trim(lst%MpiMode) /= 'NoMpi') then
288
289 ! Gathering with respect to Longitude
290 allocate(lst%allLonBeg(mmpi_npex))
291 call rpn_comm_allgather(lst%myLonBeg ,1,"mpi_integer", &
292 lst%allLonBeg,1,"mpi_integer","EW",ier)
293 if (mmpi_myid == 0) write(*,*) 'AllLonBeg =', lst%allLonBeg(:)
294
295 allocate(lst%allLonEnd(mmpi_npex))
296 call rpn_comm_allgather(lst%myLonEnd ,1,"mpi_integer", &
297 lst%allLonEnd,1,"mpi_integer","EW",ier)
298 if (mmpi_myid == 0) write(*,*) 'AllLonEnd =', lst%allLonEnd(:)
299
300 allocate(lst%allLonPerPE(mmpi_npex))
301 call rpn_comm_allgather(lst%lonPerPE ,1,"mpi_integer", &
302 lst%allLonPerPE,1,"mpi_integer","EW",ier)
303 if (mmpi_myid == 0) write(*,*) 'AllLonPerPE =', lst%allLonPerPE(:)
304
305 ! Gathering with respect to Latitude
306 allocate(lst%allLatBeg(mmpi_npey))
307 call rpn_comm_allgather(lst%myLatBeg ,1,"mpi_integer", &
308 lst%allLatBeg,1,"mpi_integer","NS",ier)
309 if (mmpi_myid == 0) write(*,*) 'AllLatBeg =', lst%allLatBeg(:)
310
311 allocate(lst%allLatEnd(mmpi_npey))
312 call rpn_comm_allgather(lst%myLatEnd ,1,"mpi_integer", &
313 lst%allLatEnd,1,"mpi_integer","NS",ier)
314 if (mmpi_myid == 0) write(*,*) 'AllLatEnd =', lst%allLatEnd(:)
315
316 allocate(lst%allLatPerPE(mmpi_npey))
317 call rpn_comm_allgather(lst%latPerPE ,1,"mpi_integer", &
318 lst%allLatPerPE,1,"mpi_integer","NS",ier)
319 if (mmpi_myid == 0) write(*,*) 'AllLatPerPE =', lst%allLatPerPE(:)
320
321 ! Gathering with respect to M
322 allocate(lst%allmBeg(mmpi_npey))
323 call rpn_comm_allgather(lst%mymBeg ,1,"mpi_integer", &
324 lst%allmBeg,1,"mpi_integer","NS",ier)
325 if (mmpi_myid == 0) write(*,*) 'AllmBeg =', lst%allmBeg(:)
326
327 allocate(lst%allmEnd(mmpi_npey))
328 call rpn_comm_allgather(lst%mymEnd ,1,"mpi_integer", &
329 lst%allmEnd,1,"mpi_integer","NS",ier)
330 if (mmpi_myid == 0) write(*,*) 'allmEnd =', lst%allmEnd(:)
331
332 allocate(lst%allmSkip(mmpi_npey))
333 call rpn_comm_allgather(lst%mymSkip ,1,"mpi_integer", &
334 lst%allmSkip,1,"mpi_integer","NS",ier)
335 if (mmpi_myid == 0) write(*,*) 'allmSkip = ', lst%allmSkip(:)
336
337 allocate(lst%allnBeg(mmpi_npex))
338 call rpn_comm_allgather(lst%mynBeg ,1,"mpi_integer", &
339 lst%allnBeg,1,"mpi_integer","EW",ier)
340 if (mmpi_myid == 0) write(*,*) 'AllnBeg =', lst%allnBeg(:)
341
342 allocate(lst%allnEnd(mmpi_npex))
343 call rpn_comm_allgather(lst%mynEnd ,1,"mpi_integer", &
344 lst%allnEnd,1,"mpi_integer","EW",ier)
345 if (mmpi_myid == 0) write(*,*) 'AllnEnd =', lst%allnEnd(:)
346
347 allocate(lst%allnSkip(mmpi_npex))
348 call rpn_comm_allgather(lst%mynSkip ,1,"mpi_integer", &
349 lst%allnSkip,1,"mpi_integer","EW",ier)
350 if (mmpi_myid == 0) write(*,*) 'AllnSkip = ', lst%allnSkip(:)
351
352 ! Gathering with respect to levels
353 call rpn_comm_allreduce(lst%myLevCount,lst%maxLevCount, &
354 1,"MPI_INTEGER","MPI_MAX","GRID",ier)
355 if (mmpi_myid == 0) write(*,*) 'MaxLevCount =',lst%maxLevCount
356
357 allocate(lst%allLevBeg(mmpi_npex))
358 call rpn_comm_allgather(lst%myLevBeg ,1,"mpi_integer", &
359 lst%allLevBeg,1,"mpi_integer","EW",ier)
360 if (mmpi_myid == 0) write(*,*) 'AllLevBeg =', lst%allLevBeg(:)
361
362 allocate(lst%allLevEnd(mmpi_npex))
363 call rpn_comm_allgather(lst%myLevEnd ,1,"mpi_integer", &
364 lst%allLevEnd,1,"mpi_integer","EW",ier)
365 if (mmpi_myid == 0) write(*,*) 'AllLevEnd =', lst%allLevEnd(:)
366
367 ! Setup mpi derived types used in transposes (only used when grid is divisible)
368 ! ... mpi_type_vector(count, blocklength, stride, ...)
369 ! ... mpi_type_create_resized(oldtype, lowerbound, extent(in bytes), newtype, ierr)
370
371 call mpi_type_size(MPI_REAL8, realSize, ierr)
372 lowerBound = 0
373
374 ! create the send type for LevToLon
375 extent = lst%maxLevCount * lst%lonPerPE * realSize
376 call mpi_type_vector(lst%latPerPE, lst%maxLevCount * lst%lonPerPE, &
377 lst%maxLevCount * lst%ni, MPI_REAL8, sendtype, ierr)
378 call mpi_type_create_resized(sendtype, lowerBound , extent, lst%sendType_LevToLon, ierr);
379 call mpi_type_commit(lst%sendType_LevToLon,ierr)
380
381 ! create the receive type for LevToLon
382 extent = lst%maxLevCount * realSize
383 call mpi_type_vector(lst%lonPerPE * lst%latPerPE , lst%maxLevCount, &
384 maxlevels_opt, MPI_REAL8, recvtype, ierr);
385 call mpi_type_create_resized(recvtype, lowerBound, extent, lst%recvType_LevToLon, ierr);
386 call mpi_type_commit(lst%recvType_LevToLon, ierr)
387
388 ! create the send type for LonToLev
389 extent = lst%maxLevCount * realSize
390 call mpi_type_vector(lst%lonPerPE * lst%latPerPE , lst%maxLevCount, &
391 maxlevels_opt, MPI_REAL8, sendtype, ierr);
392 call mpi_type_create_resized(sendtype, lowerBound, extent, lst%sendType_LonToLev, ierr);
393 call mpi_type_commit(lst%sendType_LonToLev, ierr)
394
395 ! create the recv type for LonToLev
396 extent = lst%maxLevCount * lst%lonPerPE * realSize
397 call mpi_type_vector(lst%latPerPE, lst%maxLevCount * lst%lonPerPE, &
398 lst%maxLevCount * lst%ni, MPI_REAL8, recvtype, ierr)
399 call mpi_type_create_resized(recvtype, lowerBound , extent, lst%recvType_LonToLev, ierr);
400 call mpi_type_commit(lst%recvType_LonToLev,ierr)
401
402 end if
403
404 !- 1.4 Compute the Total Wavenumber associated with weach m,n pairs and
405 ! the number of spectral element in the VAR array (nla)
406 ! FOR THE LOCAL PROCESSOR
407 allocate(Kr8fromMN(0:lst%mmax,0:lst%nmax))
408 Kr8FromMN(:,:) = -1.d0
409
410 allocate(KfromMN(0:lst%mmax,0:lst%nmax))
411 KFromMN(:,:) = -1
412
413 ! Denis et al., MWR, 2002
414 !mref = lst%mmax
415 !nref = lst%nmax
416
417 ! old var code (L. Fillion)
418 mref = lst%ni-1
419 nref = lst%nj-1
420
421 select case(trim(kreftype))
422 case ('MAX','max')
423 ! old var code (L. Fillion)
424 write(*,*)
425 write(*,*) 'lst_Setup: Using legacy (old var code) total wavenumber normalization'
426 kref = max(mref,nref)
427 case ('MIN','min')
428 ! Denis et al., MWR, 2002
429 write(*,*)
430 write(*,*) 'lst_Setup: Using Denis et al. total wavenumber normalization'
431 kref = min(mref,nref)
432 case default
433 write(*,*)
434 write(*,*) 'Unknown KREFTYPE in lst_setup : ', trim(kreftype)
435 call utl_abort('lst_setup')
436 end select
437
438 kMax = nint(lst_totalWaveNumber(lst%mmax,lst%nmax,mref,nref,kref))
439
440 if (ktrunc_in == -1) then ! no truncation case
441 lst%ktrunc = kMax
442 else
443 if (ktrunc_in > 0) then
444 lst%ktrunc = ktrunc_in
445 else
446 call utl_abort('lst_setup: invalid truncation')
447 end if
448 end if
449
450 if (lst%ktrunc >= kMax) then
451 write(*,*)
452 write(*,*) 'lst_Setup: Warning: Truncation is larger than kMax'
453 write(*,*) ' NO TRUNCATION will be applied'
454 else if (lst%ktrunc > min(lst%mmax,lst%nmax)) then
455 write(*,*)
456 if (lst%ktrunc > lst%mmax) then
457 write(*,*) 'lst_Setup: Warning: Truncation is larger than mmax only'
458 write(*,*) ' TRUNCATION will be applied only above nmax'
459 write(*,'(A,f8.1)') ' i.e., for wavelenght (in km) in y-axis smaller than ',&
460 (lst%nj*ec_ra*dlon_in/1000.0)/lst%ktrunc
461 else
462 write(*,*) 'lst_Setup: Warning: Truncation is larger than nmax only'
463 write(*,*) ' TRUNCATION will be applied only above mmax'
464 write(*,'(A,f8.1)') ' i.e., for wavelenght (in km) in x-axis smaller than ',&
465 (lst%ni*ec_ra*dlon_in/1000.0)/lst%ktrunc
466 end if
467 else
468 write(*,*)
469 write(*,*) 'lst_Setup: TRUNCATION will be applied above k = ',lst%ktrunc
470 write(*,'(A,f8.1)') ' i.e., for wavelenght (in km) in x-axis smaller than ',&
471 (lst%ni*ec_ra*dlon_in/1000.0)/lst%ktrunc
472 write(*,'(A,f8.1)') ' i.e., for wavelenght (in km) in y-axis smaller than ',&
473 (lst%nj*ec_ra*dlon_in/1000.0)/lst%ktrunc
474 end if
475
476 ila = 0
477 do n = lst%mynBeg, lst%mynEnd, lst%mynSkip
478 do m = lst%mymBeg, lst%mymEnd, lst%mymSkip
479 r = lst_totalWaveNumber(m,n,mref,nref,kref)
480 k = nint(r) ! or ceiling(r) as in Denis et al. ?
481 if (k <= lst%ktrunc) then
482 ila = ila +1
483 KFromMN(m,n) = k
484 Kr8FromMN(m,n) = r
485 end if
486 end do
487 end do
488
489 if (ila == 0) then
490 write(*,*)
491 write(*,*) 'There are no spectral elements associated to this mpi task!'
492 end if
493
494 lst%nla = ila ! Number of spectral element per phase in the VAR array
495 if (trim(lst%MpiMode) /= 'NoMpi') then
496 call rpn_comm_allreduce(lst%nla,lst%maxnla, &
497 1,"MPI_INTEGER","MPI_MAX","GRID",ier)
498 if (mmpi_myid == 0) write(*,*) 'MaxNLA =',lst%maxnla
499 end if
500
501 allocate(lst%KfromMNglb(0:lst%mmax,0:lst%nmax))
502 allocate(my_KfromMNglb(0:lst%mmax,0:lst%nmax))
503 my_KfromMNglb = 0
504 my_KFromMNglb(:,:) = KFromMN(:,:)
505 if (trim(lst%MpiMode) /= 'NoMpi') then
506 call rpn_comm_allreduce(my_KFromMNglb,lst%KFromMNglb, &
507 (lst%mmax+1)*(lst%nmax+1),"MPI_INTEGER","MPI_MAX","GRID",ier)
508 end if
509 deallocate(my_KfromMNglb)
510
511 lst%mymActiveCount=0
512 do m = lst%mymBeg, lst%mymEnd, lst%mymSkip
513 if (KfromMN(m,0) /= -1) lst%mymActiveCount = lst%mymActiveCount + 1
514 end do
515 if (trim(lst%MpiMode) /= 'NoMpi') then
516 call rpn_comm_allreduce(lst%mymActiveCount,lst%maxmActiveCount, &
517 1,"MPI_INTEGER","MPI_MAX","GRID",ier)
518 if (mmpi_myid == 0) write(*,*) 'MaxmActiveCount =',lst%maxmActiveCount
519 end if
520
521 !- 1.5 VAR spectral element ordering &
522 ! Total Wavenumbers and Weights associated with each spectral element
523 ! FOR THE LOCAL PROCESSOR
524 allocate(lst%nla_Index(0:lst%mmax,0:lst%nmax))
525
526 allocate(lst%k_r8(1:lst%nla))
527 allocate(lst%k(1:lst%nla))
528 allocate(lst%m(1:lst%nla))
529 allocate(lst%n(1:lst%nla))
530 allocate(lst%Weight(1:lst%nla))
531 allocate(lst%nePerK(0:lst%ktrunc))
532 allocate(lst%nePerKglobal(0:lst%ktrunc))
533 allocate(lst%ilaFromEK(1:lst%nla,0:lst%ktrunc))
534 allocate(lst%ilaGlobal(1:lst%nla))
535
536 lst%nla_Index(:,:) = -1
537 lst%ilaFromEK(:,:) = -1
538 lst%NEPerK(:) = 0
539 lst%NEPerKglobal(:)= 0
540
541 ila = 0
542 ilaglb = 0
543 do n = 0, lst%nmax
544 do m = 0, lst%mmax
545 k = KfromMN(m,n)
546
547 if (lst%KfromMNglb(m,n) /= -1) ilaglb = ilaglb + 1 ! Global Index
548
549 if (k /= -1) then
550 ila = ila+1
551
552 ! Internal index
553 lst%nla_Index(m,n) = ila
554
555 ! Outgoing (public) variables
556 lst%nePerK(k) = lst%nePerK(k) + 1
557 lst%ilaFromEK(lst%nePerK(k),k) = ila
558 lst%k_r8(ila) = Kr8fromMN(m,n)
559 lst%k(ila) = k
560 lst%m(ila) = m
561 lst%n(ila) = n
562 lst%ilaGlobal(ila) = ilaglb
563
564 ! Spectral coefficient weight associated with this index
565 if (m == 0 .and. n == 0) then
566 lst%Weight(ila) = 1.0d0
567 else if (m /= 0 .and. n /= 0) then
568 lst%Weight(ila) = 4.0d0
569 else
570 lst%Weight(ila) = 2.0d0
571 end if
572
573 end if
574
575 end do
576 end do
577
578 lst%nlaGlobal = ilaglb ! Number of spectral element per phase in the VAR mpi global array
579
580 if (trim(lst%MpiMode) /= 'NoMpi') then
581 call rpn_comm_allreduce(lst%nePerK, lst%nePerKglobal, lst%ktrunc+1, &
582 "mpi_integer", "mpi_sum", "GRID", ierr)
583 end if
584
585 deallocate(Kr8fromMN)
586 deallocate(KfromMN)
587
588 !- 1.7 Gridded data ordering (input/output)
589 if (present(gridDataOrder_opt)) then
590 lst%gridDataOrder = trim(gridDataOrder_opt)
591 else
592 lst%gridDataOrder = 'ijk' ! default value
593 end if
594
595 select case (trim(lst%gridDataOrder))
596 case ('ijk')
597 write(*,*) 'lst_setup: gridded data ordering = IJK'
598 case ('kij')
599 write(*,*) 'lst_setup: gridded data ordering = KIJ'
600 case default
601 write(*,*)
602 write(*,*) 'Error: gridDataOrder Unknown ', trim(gridDataOrder_opt)
603 call utl_abort('lst_setup')
604 end select
605
606 !
607 !- 2. Set factors for parseval identity
608 !
609 allocate(lst%NormFactor (lst%nla,lst%nphase))
610 allocate(lst%NormFactorAd(lst%nla,lst%nphase))
611
612 Normfactor1 = 1.0d0
613 Normfactor2 = 0.5d0 * sqrt(2.0d0)
614 Normfactor3 = 0.5d0
615 NormfactorAd1 = 1.0d0 * real((lst%ni * lst%nj),8)
616 NormfactorAd2 = sqrt(2.0d0) * real((lst%ni * lst%nj),8)
617 NormfactorAd3 = 2.0d0 * real((lst%ni * lst%nj),8)
618
619 do ila = 1,lst%nla
620
621 m = lst%m(ila)
622 n = lst%n(ila)
623
624 do p = 1, lst%nphase
625 if (p == 1) then
626 i = 2*m+1
627 j = 2*n+1
628 else if (p == 2) then
629 i = 2*m+1
630 j = 2*n+2
631 else if (p == 3) then
632 i = 2*m+2
633 j = 2*n+1
634 else if (p == 4) then
635 i = 2*m+2
636 j = 2*n+2
637 else
638 call utl_abort('lst_Setup: Error in NormFactor')
639 end if
640
641 if (i == 1 .or. j == 1) then
642 if (i == 1 .and. j == 1) then
643 factor = Normfactor1
644 factorAd = NormfactorAd1
645 else
646 factor = Normfactor2
647 factorAd = NormfactorAd2
648 end if
649 else
650 factor = Normfactor3
651 factorAd = NormfactorAd3
652 end if
653
654 lst%NormFactor (ila,p) = factor
655 lst%NormFactorAd(ila,p) = factorAd
656 end do
657
658 end do
659
660 !
661 !- 3. Set variables needed by Forward and Inverse Laplacian
662 !
663 allocate(lst%lapxy (lst%nla))
664 allocate(lst%ilapxy(lst%nla))
665
666 dlon = dlon_in
667 dx2 = (ec_ra*dlon)**2
668 fac = 2.d0/dx2
669
670 do ila = 1,lst%nla
671 ca = 2.d0*MPC_PI_R8 * lst%m(ila)
672 cp = cos(ca/lst%ni)
673 cb = 2.d0*MPC_PI_R8 * lst%n(ila)
674 cq = cos(cb/lst%nj)
675
676 lst%lapxy(ila) = fac * (cp + cq - 2.d0)
677 if (lst%lapxy(ila) /= 0.d0) then
678 lst%ilapxy(ila) = 1.d0 / lst%lapxy(ila)
679 else
680 lst%ilapxy(ila) = 0.d0
681 end if
682 end do
683
684 !
685 !- 4. Finalized
686 !
687 lst%allocated = .true.
688
689 end subroutine lst_Setup
690
691 !--------------------------------------------------------------------------
692 ! lst_totalWaveNumber
693 !--------------------------------------------------------------------------
694 function lst_totalWaveNumber(m,n,mref,nref,kref) result(r)
695 implicit none
696
697 ! Arguments:
698 integer, intent(in) :: m
699 integer, intent(in) :: n
700 integer, intent(in) :: mref
701 integer, intent(in) :: nref
702 integer, intent(in) :: kref
703 ! Result:
704 real(8) :: r
705
706 ! Locals:
707 real(8) :: a, b
708
709 a = real(m,8)/real(mref,8)
710 b = real(n,8)/real(nref,8)
711 r = real(kref,8) * sqrt((a**2) + (b**2)) ! Ellipse Shape if nref /= mref
712
713 end function lst_totalWaveNumber
714
715 !--------------------------------------------------------------------------
716 ! lst_VARTRANSFORM
717 !--------------------------------------------------------------------------
718 subroutine lst_VarTransform(lst, SpectralStateVar, GridState, &
719 TransformDirection, nk)
720 implicit none
721
722 ! Arguments:
723 type(struct_lst), intent(in) :: lst
724 integer, intent(in) :: nk ! Grid point data dimensions
725 character(len=*), intent(in) :: TransformDirection ! SpectralToGridPoint or GridPointToSpectral
726 real(8), intent(inout) :: GridState(:,:,:) ! 3D field in grid point space
727 real(8), intent(inout) :: SpectralStateVar(:,:,:) ! 3D spectral coefficients
728
729 !
730 !- 0. Some tests...
731 !
732 if (verbose) write(*,*) 'Entering lst_varTransform'
733
734 !
735 !- 1. Call the appropriate transform
736 !
737 if (trim(lst%gridDataOrder) == 'ijk') then
738 call lst_VarTransform_ijk(lst, SpectralStateVar, GridState, &
739 TransformDirection, nk)
740 else
741 call lst_VarTransform_kij(lst, SpectralStateVar, GridState, &
742 TransformDirection, nk)
743 end if
744
745 end subroutine lst_VarTransform
746
747 !--------------------------------------------------------------------------
748 ! lst_VARTRANSFORM_IJK
749 !--------------------------------------------------------------------------
750 subroutine lst_VarTransform_ijk(lst, SpectralStateVar, GridState, &
751 TransformDirection, nk)
752 implicit none
753
754 ! Arguments:
755 type(struct_lst), intent(in) :: lst
756 integer, intent(in) :: nk ! Grid point data dimensions
757 character(len=*), intent(in) :: TransformDirection ! SpectralToGridPoint or GridPointToSpectral
758 real(8), intent(inout) :: GridState(lst%myLonBeg:lst%myLonEnd,lst%myLatBeg:lst%myLatEnd,nk) ! 3D field in grid point space
759 real(8), intent(inout) :: SpectralStateVar(:,:,:) ! 3D spectral coefficients
760
761 ! Locals:
762 integer :: ni_l, nj_l, nip_l, njp_l
763 integer :: iStart, iEnd, jStart, jEnd, kStart, kEnd
764 real(8), allocatable :: Step0(:,:,:)
765 real(8), allocatable :: Step1(:,:,:)
766 real(8), allocatable :: Step2(:,:,:)
767 real(8), allocatable :: Step3(:,:,:)
768 character(len=1) :: TransformAxe
769
770 !
771 !- 1. Data pre-processing (Input -> Step0)
772 !
773 if (verbose) write(*,*) 'Entering lst_varTransform_ijk'
774
775 !- 1.1 Settings and Data Selection
776
777 if (trim(TransformDirection) == 'GridPointToSpectral') then
778 iStart = 1
779 iEnd = lst%ni
780 jStart = lst%myLatBeg
781 jEnd = lst%myLatEnd
782 if (trim(lst%MpiMode) == 'NoMpi') then
783 kStart= 1
784 kEnd = nk
785 else
786 kStart= lst%myLevBeg
787 kEnd = lst%myLevEnd
788 end if
789 allocate(Step0(iStart:iEnd,jStart:jEnd,kStart:kEnd))
790 if (trim(lst%MpiMode) == 'NoMpi') then
791 Step0(:,:,:) = GridState(:,:,:)
792 else
793 call transpose2d_LonToLev(Step0, & ! OUT
794 GridState, nk, lst) ! IN
795 end if
796 end if
797
798 !
799 !- 1. First pass (Step0 -> Step1)
800 !
801
802 !- 1.1 Settings and Data Selection
803 select case (trim(TransformDirection))
804 case ('GridPointToSpectral')
805 TransformAxe = 'i'
806 ni_l = lst%ni
807 nip_l = lst%nip
808 nj_l = lst%latPerPE
809 njp_l = 0
810 if (trim(lst%MpiMode) == 'NoMpi') then
811 kStart= 1
812 kEnd = nk
813 else
814 kStart= lst%myLevBeg
815 kEnd = lst%myLevEnd
816 end if
817 case ('SpectralToGridPoint')
818 TransformAxe = 'j'
819 ni_l = 2*lst%mymCount
820 nip_l = 0
821 nj_l = lst%nj
822 njp_l = lst%njp
823 if (trim(lst%MpiMode) == 'NoMpi') then
824 kStart= 1
825 kEnd = nk
826 else
827 kStart= lst%myLevBeg
828 kEnd = lst%myLevEnd
829 end if
830 case default
831 write(*,*)
832 write(*,*) 'Error: TranformDirection Unknown ', trim(TransformDirection)
833 call utl_abort('lst_VarTransform')
834 end select
835
836 allocate(Step1(ni_l+nip_l,nj_l+njp_l,kStart:kEnd))
837
838 !- 1.2 Spectral transform
839 if (trim(TransformDirection) == 'SpectralToGridPoint') then
840 if (trim(lst%MpiMode) == 'NoMpi') then
841 call lst_ReshapeTrunc(Step1, & ! OUT
842 SpectralStateVar, & ! IN
843 'ToRPN', kStart, kEnd, lst) ! IN
844 else
845 call transpose2d_NToLev(Step1, & ! OUT
846 SpectralStateVar, nk, lst) ! IN
847 end if
848 else
849 Step1(1:lst%ni,1:lst%latPerPE,:) = Step0(1:lst%ni,lst%myLatBeg:lst%myLatEnd,:)
850 deallocate(Step0)
851 end if
852
853 call lst_transform1d(Step1, & ! INOUT
854 TransformDirection, & ! IN
855 TransformAxe, & ! IN
856 ni_l, nj_l, nip_l, njp_l, & ! IN
857 kStart, kEnd) ! IN
858
859 !
860 !- 2.0 Second pass (Step1 -> Step2)
861 !
862
863 !- 2.1 Settings
864 if (trim(TransformDirection) == 'GridPointToSpectral') then
865 TransformAxe = 'j'
866 ni_l = 2*lst%mymCount
867 nip_l = 0
868 nj_l = lst%nj
869 njp_l = lst%njp
870 if (trim(lst%MpiMode) == 'NoMpi') then
871 kStart= 1
872 kEnd = nk
873 else
874 kStart= lst%myLevBeg
875 kEnd = lst%myLevEnd
876 end if
877 else
878 TransformAxe = 'i'
879 ni_l = lst%ni
880 nip_l = lst%nip
881 nj_l = lst%latPerPE
882 njp_l = 0
883 if (trim(lst%MpiMode) == 'NoMpi') then
884 kStart= 1
885 kEnd = nk
886 else
887 kStart= lst%myLevBeg
888 kEnd = lst%myLevEnd
889 end if
890 end if
891
892 allocate(Step2(ni_l+nip_l,nj_l+njp_l,kStart:kEnd))
893
894 !- 2.2 Communication between processors
895
896 if (trim(TransformDirection) == 'GridPointToSpectral') then
897 if (trim(lst%MpiMode) == 'NoMpi') then
898 Step2(:,1:lst%nj,:) = Step1(:,1:lst%nj,:)
899 else
900 call transpose2d_LatToM(Step2, & ! OUT
901 Step1, lst) ! IN
902 end if
903 else
904 if (trim(lst%MpiMode) == 'NoMpi') then
905 Step2(:,1:lst%nj,:) = Step1(:,1:lst%nj,:)
906 else
907 call transpose2d_MToLat(Step2, & ! OUT
908 Step1, lst) ! IN
909 end if
910 end if
911
912 deallocate(Step1)
913
914 !- 2.3 Spectral Transform
915 call lst_transform1d(Step2, & ! INOUT
916 TransformDirection, & ! IN
917 TransformAxe, & ! IN
918 ni_l, nj_l, nip_l, njp_l, & ! IN
919 kStart, kEnd) ! IN
920
921 !
922 !- 3.0 Post-processing (Step2 -> Step3 -> Output)
923 !
924
925 select case (trim(TransformDirection))
926 case ('GridPointToSpectral')
927 iStart = 1
928 iEnd = 2*lst%mymCount
929 jStart = 1
930 jEnd = 2*lst%mynCount
931 kStart= 1
932 kEnd = nk
933 case ('SpectralToGridPoint')
934 iStart = 1
935 iEnd = lst%lonPerPE
936 jStart = 1
937 jEnd = lst%latPerPE
938 kStart = 1
939 kEnd = nk
940 case default
941 write(*,*)
942 write(*,*) 'Error: TranformDirection Unknown ', trim(TransformDirection)
943 call utl_abort('lst_VarTransform')
944 end select
945
946 if (trim(TransformDirection) == 'GridPointToSpectral') then
947
948 ! Communication between processors (Truncation (if applicable) will occur in this step)
949 if (trim(lst%MpiMode) == 'NoMpi') then
950 call lst_ReshapeTrunc(Step2, & ! IN
951 SpectralStateVar, & ! OUT
952 'ToVAR', kStart, kEnd, lst) ! IN
953 else
954 call transpose2d_LevToN(SpectralStateVar, & ! OUT
955 Step2, nk, lst) ! IN
956 end if
957
958 else
959
960 allocate(Step3(iStart:iEnd,jStart:jEnd,kStart:kEnd))
961
962 ! Communication between processors
963 if (trim(lst%MpiMode) == 'NoMpi') then
964 Step3(:,:,:) = Step2(1:lst%ni,:,:)
965 else
966 call transpose2d_LevToLon(Step3, & ! OUT
967 Step2(1:lst%ni,:,:), nk, lst) ! IN
968 end if
969
970 GridState(lst%myLonBeg:lst%myLonEnd,lst%myLatBeg:lst%myLatEnd,:) = Step3(1:lst%lonPerPE,1:lst%latPerPE,:)
971
972 deallocate(Step3)
973
974 end if
975
976 deallocate(Step2)
977
978 end subroutine lst_VarTransform_ijk
979
980 !--------------------------------------------------------------------------
981 ! lst_VARTRANSFORM_KIJ
982 !--------------------------------------------------------------------------
983 subroutine lst_VarTransform_kij(lst, SpectralStateVar, GridState, &
984 TransformDirection, nk)
985 implicit none
986
987 ! Arguments:
988 type(struct_lst), intent(in) :: lst
989 integer, intent(in) :: nk ! Grid point data dimensions
990 character(len=*), intent(in) :: TransformDirection ! SpectralToGridPoint or GridPointToSpectral
991 real(8), intent(inout) :: GridState(nk,lst%myLonBeg:lst%myLonEnd,lst%myLatBeg:lst%myLatEnd) ! 3D field in grid point space
992 real(8), intent(inout) :: SpectralStateVar(:,:,:) ! 3D spectral coefficients
993
994 ! Locals:
995 integer :: ni_l, nj_l, nip_l, njp_l
996 integer :: iStart, iEnd, jStart, jEnd, kStart, kEnd
997 real(8), allocatable :: Step0(:,:,:)
998 real(8), allocatable :: Step1(:,:,:)
999 real(8), allocatable :: Step2(:,:,:)
1000 real(8), allocatable :: Step3(:,:,:)
1001 character(len=1) :: TransformAxe
1002
1003 !
1004 !- 1. Data pre-processing (Input -> Step0)
1005 !
1006 if (verbose) write(*,*) 'Entering lst_varTransform_kij'
1007
1008 !- 1.1 Settings and Data Selection
1009
1010 if (trim(TransformDirection) == 'GridPointToSpectral') then
1011 iStart = 1
1012 iEnd = lst%ni
1013 jStart = lst%myLatBeg
1014 jEnd = lst%myLatEnd
1015 if (trim(lst%MpiMode) == 'NoMpi') then
1016 kStart= 1
1017 kEnd = nk
1018 else
1019 kStart= lst%myLevBeg
1020 kEnd = lst%myLevEnd
1021 end if
1022 allocate(Step0(kStart:kEnd,iStart:iEnd,jStart:jEnd))
1023 if (trim(lst%MpiMode) == 'NoMpi') then
1024 Step0(:,:,:) = GridState(:,:,:)
1025 else
1026 if(lst%lonLatDivisible) then
1027 call transpose2d_LonToLev_kij_mpitypes(Step0, & ! OUT
1028 GridState, nk, lst) ! IN
1029 else
1030 call transpose2d_LonToLev_kij(Step0, & ! OUT
1031 GridState, nk, lst) ! IN
1032 end if
1033 end if
1034 end if
1035
1036 !
1037 !- 1. First pass (Step0 -> Step1)
1038 !
1039
1040 !- 1.1 Settings and Data Selection
1041 select case (trim(TransformDirection))
1042 case ('GridPointToSpectral')
1043 TransformAxe = 'i'
1044 ni_l = lst%ni
1045 nip_l = lst%nip
1046 nj_l = lst%latPerPE
1047 njp_l = 0
1048 if (trim(lst%MpiMode) == 'NoMpi') then
1049 kStart= 1
1050 kEnd = nk
1051 else
1052 kStart= lst%myLevBeg
1053 kEnd = lst%myLevEnd
1054 end if
1055 case ('SpectralToGridPoint')
1056 TransformAxe = 'j'
1057 ni_l = 2*lst%mymCount
1058 nip_l = 0
1059 nj_l = lst%nj
1060 njp_l = lst%njp
1061 if (trim(lst%MpiMode) == 'NoMpi') then
1062 kStart= 1
1063 kEnd = nk
1064 else
1065 kStart= lst%myLevBeg
1066 kEnd = lst%myLevEnd
1067 end if
1068 case default
1069 write(*,*)
1070 write(*,*) 'Error: TranformDirection Unknown ', trim(TransformDirection)
1071 call utl_abort('lst_VarTransform')
1072 end select
1073
1074 allocate(Step1(kStart:kEnd,ni_l+nip_l,nj_l+njp_l))
1075
1076 !- 1.2 Spectral transform
1077 if (trim(TransformDirection) == 'SpectralToGridPoint') then
1078 if (trim(lst%MpiMode) == 'NoMpi') then
1079 call lst_ReshapeTrunc_kij(Step1, & ! OUT
1080 SpectralStateVar, & ! IN
1081 'ToRPN', kStart, kEnd, lst) ! IN
1082 else
1083 call transpose2d_NToLev_kij(Step1, & ! OUT
1084 SpectralStateVar, nk, lst) ! IN
1085 end if
1086 else
1087 Step1(:,1:lst%ni,1:lst%latPerPE) = Step0(:,1:lst%ni,lst%myLatBeg:lst%myLatEnd)
1088 deallocate(Step0)
1089 end if
1090
1091 call lst_transform1d_kij(Step1, & ! INOUT
1092 TransformDirection, & ! IN
1093 TransformAxe, & ! IN
1094 ni_l, nj_l, nip_l, njp_l, & ! IN
1095 kStart, kEnd) ! IN
1096
1097 !
1098 !- 2.0 Second pass (Step1 -> Step2)
1099 !
1100
1101 !- 2.1 Settings
1102 if (trim(TransformDirection) == 'GridPointToSpectral') then
1103 TransformAxe = 'j'
1104 ni_l = 2*lst%mymCount
1105 nip_l = 0
1106 nj_l = lst%nj
1107 njp_l = lst%njp
1108 if (trim(lst%MpiMode) == 'NoMpi') then
1109 kStart= 1
1110 kEnd = nk
1111 else
1112 kStart= lst%myLevBeg
1113 kEnd = lst%myLevEnd
1114 end if
1115 else
1116 TransformAxe = 'i'
1117 ni_l = lst%ni
1118 nip_l = lst%nip
1119 nj_l = lst%latPerPE
1120 njp_l = 0
1121 if (trim(lst%MpiMode) == 'NoMpi') then
1122 kStart= 1
1123 kEnd = nk
1124 else
1125 kStart= lst%myLevBeg
1126 kEnd = lst%myLevEnd
1127 end if
1128 end if
1129
1130 allocate(Step2(kStart:kEnd,ni_l+nip_l,nj_l+njp_l))
1131
1132 !- 2.2 Communication between processors
1133
1134 if (trim(TransformDirection) == 'GridPointToSpectral') then
1135 if (trim(lst%MpiMode) == 'NoMpi') then
1136 Step2(:,1:lst%nj,:) = Step1(:,1:lst%nj,:)
1137 else
1138 call transpose2d_LatToM_kij(Step2, & ! OUT
1139 Step1, lst) ! IN
1140 end if
1141 else
1142 if (trim(lst%MpiMode) == 'NoMpi') then
1143 Step2(:,1:lst%nj,:) = Step1(:,1:lst%nj,:)
1144 else
1145 call transpose2d_MToLat_kij(Step2, & ! OUT
1146 Step1, lst) ! IN
1147 end if
1148 end if
1149
1150 deallocate(Step1)
1151
1152 !- 2.3 Spectral Transform
1153 call lst_transform1d_kij(Step2, & ! INOUT
1154 TransformDirection, & ! IN
1155 TransformAxe, & ! IN
1156 ni_l, nj_l, nip_l, njp_l, & ! IN
1157 kStart, kEnd) ! IN
1158
1159 !
1160 !- 3.0 Post-processing (Step2 -> Step3 -> Output)
1161 !
1162
1163 select case (trim(TransformDirection))
1164 case ('GridPointToSpectral')
1165 iStart = 1
1166 iEnd = 2*lst%mymCount
1167 jStart = 1
1168 jEnd = 2*lst%mynCount
1169 kStart= 1
1170 kEnd = nk
1171 case ('SpectralToGridPoint')
1172 iStart = 1
1173 iEnd = lst%lonPerPE
1174 jStart = 1
1175 jEnd = lst%latPerPE
1176 kStart = 1
1177 kEnd = nk
1178 case default
1179 write(*,*)
1180 write(*,*) 'Error: TranformDirection Unknown ', trim(TransformDirection)
1181 call utl_abort('lst_VarTransform')
1182 end select
1183
1184 if (trim(TransformDirection) == 'GridPointToSpectral') then
1185
1186 ! Communication between processors (Truncation (if applicable) will occur in this step)
1187 if (trim(lst%MpiMode) == 'NoMpi') then
1188 call lst_ReshapeTrunc_kij(Step2, & ! IN
1189 SpectralStateVar, & ! OUT
1190 'ToVAR', kStart, kEnd, lst) ! IN
1191 else
1192 call transpose2d_LevToN_kij(SpectralStateVar, & ! OUT
1193 Step2, nk, lst) ! IN
1194 end if
1195
1196 else
1197
1198 allocate(Step3(kStart:kEnd,iStart:iEnd,jStart:jEnd))
1199
1200 ! Communication between processors
1201 if (trim(lst%MpiMode) == 'NoMpi') then
1202 Step3(:,:,:) = Step2(:,1:lst%ni,:)
1203 else
1204 if(lst%lonLatDivisible) then
1205 call transpose2d_LevToLon_kij_mpitypes(Step3, & ! OUT
1206 Step2(:,1:lst%ni,:), nk, lst) ! IN
1207 else
1208 call transpose2d_LevToLon_kij(Step3, & ! OUT
1209 Step2(:,1:lst%ni,:), nk, lst) ! IN
1210 end if
1211 end if
1212
1213 GridState(:,lst%myLonBeg:lst%myLonEnd,lst%myLatBeg:lst%myLatEnd) = Step3(:,1:lst%lonPerPE,1:lst%latPerPE)
1214
1215 deallocate(Step3)
1216
1217 end if
1218
1219 deallocate(Step2)
1220
1221 end subroutine lst_VarTransform_kij
1222
1223 !--------------------------------------------------------------------------
1224 ! lst_TRANSFORM1D
1225 !--------------------------------------------------------------------------
1226 subroutine lst_transform1d(Field3d, &
1227 TransformDirection, &
1228 TransformAxe, &
1229 ni_l, nj_l, nip_l, njp_l, kStart, kEnd)
1230 implicit none
1231
1232 ! Arguments:
1233 integer, intent(in) :: ni_l ! Grid point data dimensions
1234 integer, intent(in) :: nj_l ! Grid point data dimensions
1235 integer, intent(in) :: kStart ! Grid point data dimensions
1236 integer, intent(in) :: kEnd ! Grid point data dimensions
1237 integer, intent(in) :: nip_l ! Extra point in spectral space
1238 integer, intent(in) :: njp_l ! Extra point in spectral space
1239 character(len=*), intent(in) :: TransformDirection ! SpectralToGridPoint or GridPointToSpectral
1240 character(len=*), intent(in) :: TransformAxe ! 'i' or 'j'
1241 real(8), intent(inout) :: Field3d(1:ni_l+nip_l,1:nj_l+njp_l,kStart:kEnd) ! InOut 3D field
1242
1243 ! Locals:
1244 integer :: way
1245 integer :: axe, n, nlot, nfact, np, lot, nk
1246
1247 !
1248 !- 1. Set some options
1249 !
1250 if (verbose) write(*,*) 'Entering lst_transform1d'
1251 call utl_tmg_start(151,'low-level--lst_fft')
1252
1253 !- 1.1 Transform Direction
1254 select case (trim(TransformDirection))
1255 case ('GridPointToSpectral')
1256 way = -1
1257 case ('SpectralToGridPoint')
1258 way = +1
1259 case default
1260 write(*,*)
1261 write(*,*) 'Error: TranformDirection Unknown ', trim(TransformDirection)
1262 call utl_abort('lst_VarTransform')
1263 end select
1264
1265 nk = kEnd - kStart + 1
1266
1267 !
1268 !- 2. Do the transforms in one direction
1269 !
1270 select case (trim(TransformAxe))
1271 case ('i')
1272 !- 2.1 First pass --> Along INDEX "I"
1273 axe = 0
1274 n = ni_l
1275 np = nip_l
1276 nlot= nj_l
1277 case ('j')
1278 !- 2.2 Second pass --> Along INDEX "J"
1279 axe = 1
1280 n = nj_l
1281 np = njp_l
1282 nlot= ni_l
1283 case default
1284 write(*,*)
1285 write(*,*) 'Error: TranformAxe Unknown ', trim(TransformAxe)
1286 call utl_abort('lst_VarTransform')
1287 end select
1288
1289 !- 1.2 Fast or Slow Fourier Transform ?
1290 nfact = n
1291 call ngfft(nfact) ! INOUT
1292
1293 if (nfact == n) then
1294 call setfft8(n) ! IN
1295 else
1296 call utl_abort('lst_VarTransform: This module can only handle fast sin&cos FFT')
1297 end if
1298
1299 select case (trim(TransformAxe))
1300 case ('i')
1301 !$OMP PARALLEL DO PRIVATE(lot)
1302 do lot = 1, nlot
1303 call ffft8(Field3d(:,lot,:), 1, n+np, nk, way)
1304 end do
1305 !$OMP END PARALLEL DO
1306 case ('j')
1307 !$OMP PARALLEL DO PRIVATE(lot)
1308 do lot = 1, nlot
1309 call ffft8(Field3d(lot,:,:), 1, n+np, nk, way)
1310 end do
1311 !$OMP END PARALLEL DO
1312 end select
1313
1314 !* subroutine ffft8( a, inc, jump, lot, isign )
1315 !* a is the array containing input & output data
1316 !* inc is the increment within each data 'vector'
1317 !* (e.g. inc=1 for consecutively stored data)
1318 !* jump is the increment between the start of each data vector
1319 !* lot is the number of data vectors
1320 !* isign = +1 for transform from spectral to gridpoint
1321 !* = -1 for transform from gridpoint to spectral
1322
1323 call utl_tmg_stop(151)
1324
1325 end subroutine lst_transform1d
1326
1327 !--------------------------------------------------------------------------
1328 ! lst_TRANSFORM1D_KIJ
1329 !--------------------------------------------------------------------------
1330 subroutine lst_transform1d_kij(Field3d, &
1331 TransformDirection, &
1332 TransformAxe, &
1333 ni_l, nj_l, nip_l, njp_l, kStart, kEnd)
1334 implicit none
1335
1336 ! Arguments:
1337 integer, intent(in) :: ni_l ! Grid point data dimensions
1338 integer, intent(in) :: nj_l ! Grid point data dimensions
1339 integer, intent(in) :: kStart ! Grid point data dimensions
1340 integer, intent(in) :: kEnd ! Grid point data dimensions
1341 integer, intent(in) :: nip_l ! Extra point in spectral space
1342 integer, intent(in) :: njp_l ! Extra point in spectral space
1343 character(len=*), intent(in) :: TransformDirection ! SpectralToGridPoint or GridPointToSpectral
1344 character(len=*), intent(in) :: TransformAxe ! 'i' or 'j'
1345 real(8), intent(inout) :: Field3d(kStart:kEnd,1:ni_l+nip_l,1:nj_l+njp_l) ! InOut 3D field
1346
1347 ! Locals:
1348 integer :: way
1349 integer :: axe, n, nlot, nfact, np, lot, nk, k
1350
1351 !
1352 !- 1. Set some options
1353 !
1354 if (verbose) write(*,*) 'Entering lst_transform1d_kij'
1355 call utl_tmg_start(151,'low-level--lst_fft')
1356
1357 !- 1.1 Transform Direction
1358 select case (trim(TransformDirection))
1359 case ('GridPointToSpectral')
1360 way = -1
1361 case ('SpectralToGridPoint')
1362 way = +1
1363 case default
1364 write(*,*)
1365 write(*,*) 'Error: TranformDirection Unknown ', trim(TransformDirection)
1366 call utl_abort('lst_VarTransform')
1367 end select
1368
1369 nk = kEnd - kStart + 1
1370
1371 !
1372 !- 2. Do the transforms in one direction
1373 !
1374 select case (trim(TransformAxe))
1375 case ('i')
1376 !- 2.1 First pass --> Along INDEX "I"
1377 axe = 0
1378 n = ni_l
1379 np = nip_l
1380 nlot= nj_l
1381 case ('j')
1382 !- 2.2 Second pass --> Along INDEX "J"
1383 axe = 1
1384 n = nj_l
1385 np = njp_l
1386 nlot= ni_l
1387 case default
1388 write(*,*)
1389 write(*,*) 'Error: TranformAxe Unknown ', trim(TransformAxe)
1390 call utl_abort('lst_VarTransform')
1391 end select
1392
1393 !- 1.2 Fast or Slow Fourier Transform ?
1394 nfact = n
1395 call ngfft(nfact) ! INOUT
1396
1397 if (nfact == n) then
1398 call setfft8(n) ! IN
1399 else
1400 call utl_abort('lst_VarTransform: This module can only handle fast sin&cos FFT')
1401 end if
1402
1403 select case (trim(TransformAxe))
1404 case ('i')
1405 !$OMP PARALLEL DO PRIVATE(lot,k)
1406 do lot = 1, nlot
1407 call ffft8(Field3d(:,:,lot), nk, 1, nk, way)
1408 end do
1409 !$OMP END PARALLEL DO
1410 case ('j')
1411 !$OMP PARALLEL DO PRIVATE(lot,k)
1412 do lot = 1, nlot
1413 call ffft8(Field3d(:,lot,:), nk, 1, nk, way)
1414 end do
1415 !$OMP END PARALLEL DO
1416 end select
1417
1418 !* subroutine ffft8( a, inc, jump, lot, isign )
1419 !* a is the array containing input & output data
1420 !* inc is the increment within each data 'vector'
1421 !* (e.g. inc=1 for consecutively stored data)
1422 !* jump is the increment between the start of each data vector
1423 !* lot is the number of data vectors
1424 !* isign = +1 for transform from spectral to gridpoint
1425 !* = -1 for transform from gridpoint to spectral
1426
1427 call utl_tmg_stop(151)
1428
1429 end subroutine lst_transform1d_kij
1430
1431 !--------------------------------------------------------------------------
1432 ! lst_Transpose2d_LonToLev
1433 !--------------------------------------------------------------------------
1434 subroutine transpose2d_LonToLev(gd_out, gd_in, nk, lst)
1435 implicit none
1436
1437 ! Arguments:
1438 type(struct_lst), intent(in) :: lst
1439 integer, intent(in) :: nk
1440 real(8), intent(in) :: gd_in(lst%myLonBeg:lst%myLonEnd, lst%myLatBeg:lst%myLatEnd, nk)
1441 real(8), intent(out) :: gd_out(lst%ni, lst%myLatBeg:lst%myLatEnd, lst%myLevBeg:lst%myLevEnd)
1442
1443 ! Locals:
1444 real(8) :: gd_send(lst%lonPerPEmax, lst%latPerPEmax, lst%maxLevCount, mmpi_npex)
1445 real(8) :: gd_recv(lst%lonPerPEmax, lst%latPerPEmax, lst%maxLevCount, mmpi_npex)
1446 integer :: yourid, nsize, ierr, levIndex, levIndex2
1447
1448 if (verbose) write(*,*) 'Entering transpose2d_LonToLev'
1449 call rpn_comm_barrier("GRID",ierr)
1450
1451 call utl_tmg_start(155,'low-level--lst_transpose_LEVtoLON')
1452
1453 !$OMP PARALLEL DO PRIVATE(yourid,levIndex,levIndex2)
1454 do yourid = 0, (mmpi_npex-1)
1455 gd_send(:,:,:,yourid+1) = 0.0d0
1456 do levIndex = lst%allLevBeg(yourid+1), lst%allLevEnd(yourid+1)
1457 levIndex2 = levIndex-lst%allLevBeg(yourid+1)+1
1458 gd_send(1:lst%lonPerPE,1:lst%latPerPE,levIndex2,yourid+1) = gd_in(:,:,levIndex)
1459 end do
1460 end do
1461 !$OMP END PARALLEL DO
1462
1463 nsize = lst%lonPerPEmax * lst%maxLevCount * lst%latPerPEmax
1464 if (mmpi_npex > 1) then
1465 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
1466 gd_recv,nsize,"mpi_double_precision","EW",ierr)
1467 else
1468 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
1469 end if
1470
1471 !$OMP PARALLEL DO PRIVATE(yourid,levIndex,levIndex2)
1472 do levIndex = lst%myLevBeg, lst%myLevEnd
1473 levIndex2 = levIndex - lst%myLevBeg + 1
1474 do yourid = 0, (mmpi_npex-1)
1475 gd_out(lst%allLonBeg(yourid+1):lst%allLonEnd(yourid+1),:,levIndex) = &
1476 gd_recv(1:lst%allLonPerPE(yourid+1),1:lst%latPerPE,levIndex2,yourid+1)
1477 end do
1478 end do
1479 !$OMP END PARALLEL DO
1480
1481 call utl_tmg_stop(155)
1482
1483 end subroutine transpose2d_LonToLev
1484
1485 !--------------------------------------------------------------------------
1486 ! lst_Transpose2d_LonToLev_kij_mpitypes
1487 !--------------------------------------------------------------------------
1488 subroutine transpose2d_LonToLev_kij_mpitypes(gd_out, gd_in, nk, lst)
1489 implicit none
1490
1491 ! Arguments:
1492 type(struct_lst), intent(in) :: lst
1493 integer, intent(in) :: nk
1494 real(8), intent(in) :: gd_in (nk,lst%myLonBeg:lst%myLonEnd, lst%myLatBeg:lst%myLatEnd)
1495 real(8), intent(out) :: gd_out(lst%myLevBeg:lst%myLevEnd,lst%ni, lst%myLatBeg:lst%myLatEnd)
1496
1497 ! Locals:
1498 integer :: nsize, ierr
1499
1500 if (verbose) write(*,*) 'Entering transpose2d_LonToLev_kij'
1501 call rpn_comm_barrier("GRID",ierr)
1502
1503 call utl_tmg_start(155,'low-level--lst_transpose_LEVtoLON')
1504
1505 nsize = lst%lonPerPE * lst%maxLevCount * lst%latPerPE
1506 if (mmpi_npex > 1) then
1507 call mpi_alltoall(gd_in, 1, lst%sendType_LonToLev, &
1508 gd_out, 1, lst%recvType_LonToLev, mmpi_comm_EW, ierr)
1509 else
1510 gd_out(:,:,:) = gd_in(:,:,:)
1511 end if
1512
1513 call utl_tmg_stop(155)
1514
1515 end subroutine transpose2d_LonToLev_kij_mpitypes
1516
1517 !--------------------------------------------------------------------------
1518 ! lst_Transpose2d_LonToLev_kij
1519 !--------------------------------------------------------------------------
1520 subroutine transpose2d_LonToLev_kij(gd_out, gd_in, nk, lst)
1521 implicit none
1522
1523 ! Arguments:
1524 type(struct_lst), intent(in) :: lst
1525 integer, intent(in) :: nk
1526 real(8), intent(in) :: gd_in (nk,lst%myLonBeg:lst%myLonEnd, lst%myLatBeg:lst%myLatEnd)
1527 real(8), intent(out) :: gd_out(lst%myLevBeg:lst%myLevEnd,lst%ni, lst%myLatBeg:lst%myLatEnd)
1528
1529 ! Locals:
1530 real(8) :: gd_send(lst%maxLevCount,lst%lonPerPEmax, lst%latPerPEmax, mmpi_npex)
1531 real(8) :: gd_recv(lst%maxLevCount,lst%lonPerPEmax, lst%latPerPEmax, mmpi_npex)
1532 integer :: yourid, nsize, ierr, levIndex, levIndex2, latIndex, latIndex2, lonIndex, lonIndex2
1533
1534 if (verbose) write(*,*) 'Entering transpose2d_LonToLev_kij'
1535 call rpn_comm_barrier("GRID",ierr)
1536
1537 call utl_tmg_start(155,'low-level--lst_transpose_LEVtoLON')
1538
1539 !$OMP PARALLEL DO PRIVATE(yourid,latIndex,latIndex2,levIndex,levIndex2,lonIndex,lonIndex2)
1540 do yourid = 0, (mmpi_npex-1)
1541 gd_send(:,:,:,yourid+1) = 0.0d0
1542 do latIndex = lst%myLatBeg, lst%myLatEnd
1543 latIndex2 = latIndex - lst%myLatBeg + 1
1544 do lonIndex = lst%myLonBeg, lst%myLonEnd
1545 lonIndex2 = lonIndex - lst%myLonBeg + 1
1546 do levIndex = lst%allLevBeg(yourid+1), lst%allLevEnd(yourid+1)
1547 levIndex2 = levIndex-lst%allLevBeg(yourid+1)+1
1548 gd_send(levIndex2,lonIndex2,latIndex2,yourid+1) = gd_in(levIndex,lonIndex,latIndex)
1549 end do
1550 end do
1551 end do
1552 end do
1553 !$OMP END PARALLEL DO
1554
1555 nsize = lst%lonPerPEmax * lst%maxLevCount * lst%latPerPEmax
1556 if (mmpi_npex > 1) then
1557 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
1558 gd_recv,nsize,"mpi_double_precision","EW",ierr)
1559 else
1560 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
1561 end if
1562
1563 !$OMP PARALLEL DO PRIVATE(yourid,levIndex,levIndex2,lonIndex,lonIndex2,latIndex,latIndex2)
1564 do yourid = 0, (mmpi_npex-1)
1565 do latIndex = lst%myLatBeg, lst%myLatEnd
1566 latIndex2 = latIndex - lst%myLatBeg + 1
1567 do lonIndex = lst%allLonBeg(yourid+1), lst%allLonEnd(yourid+1)
1568 lonIndex2 = lonIndex - lst%allLonBeg(yourid+1) + 1
1569 do levIndex = lst%myLevBeg, lst%myLevEnd
1570 levIndex2 = levIndex - lst%myLevBeg + 1
1571 gd_out(levIndex,lonIndex,latIndex) = gd_recv(levIndex2,lonIndex2,latIndex2,yourid+1)
1572 end do
1573 end do
1574 end do
1575 end do
1576 !$OMP END PARALLEL DO
1577
1578 call utl_tmg_stop(155)
1579
1580 end subroutine transpose2d_LonToLev_kij
1581
1582 !--------------------------------------------------------------------------
1583 ! lst_Transpose2d_LevToLon
1584 !--------------------------------------------------------------------------
1585 subroutine transpose2d_LevToLon(gd_out,gd_in,nk,lst)
1586 implicit none
1587
1588 ! Arguments:
1589 type(struct_lst), intent(in) :: lst
1590 integer, intent(in) :: nk
1591 real(8), intent(out) :: gd_out(lst%myLonBeg:lst%myLonEnd, lst%myLatBeg:lst%myLatEnd, nk)
1592 real(8), intent(in) :: gd_in(lst%ni, lst%myLatBeg:lst%myLatEnd, lst%myLevBeg:lst%myLevEnd)
1593
1594 ! Locals:
1595 real(8) :: gd_send(lst%lonPerPEmax, lst%latPerPEmax, lst%maxLevCount, mmpi_npex)
1596 real(8) :: gd_recv(lst%lonPerPEmax, lst%latPerPEmax, lst%maxLevCount, mmpi_npex)
1597 integer :: yourid, nsize, ierr, levIndex, levIndex2
1598
1599 if (verbose) write(*,*) 'Entering transpose2d_LevToLon'
1600 call rpn_comm_barrier("GRID",ierr)
1601
1602 call utl_tmg_start(155,'low-level--lst_transpose_LEVtoLON')
1603
1604 !$OMP PARALLEL DO PRIVATE(yourid,levIndex,levIndex2)
1605 do levIndex = lst%myLevBeg, lst%myLevEnd
1606 levIndex2 = levIndex - lst%myLevBeg + 1
1607 do yourid = 0, (mmpi_npex-1)
1608 gd_send(:,:,levIndex2,yourid+1) = 0.0d0
1609 gd_send(1:lst%allLonPerPE(yourid+1),1:lst%latPerPE,levIndex2,yourid+1) = &
1610 gd_in(lst%allLonBeg(yourid+1):lst%allLonEnd(yourid+1),:,levIndex)
1611 end do
1612 end do
1613 !$OMP END PARALLEL DO
1614
1615 nsize = lst%lonPerPEmax * lst%maxLevCount * lst%latPerPEmax
1616 if (mmpi_npex > 1) then
1617 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
1618 gd_recv,nsize,"mpi_double_precision","EW",ierr)
1619 else
1620 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
1621 end if
1622
1623 !$OMP PARALLEL DO PRIVATE(yourid,levIndex,levIndex2)
1624 do yourid = 0, (mmpi_npex-1)
1625 do levIndex=lst%allLevBeg(yourid+1),lst%allLevEnd(yourid+1)
1626 levIndex2=levIndex-lst%allLevBeg(yourid+1)+1
1627 gd_out(:,:,levIndex) = gd_recv(1:lst%lonPerPE,1:lst%latPerPE,levIndex2,yourid+1)
1628 end do
1629 end do
1630 !$OMP END PARALLEL DO
1631
1632 call utl_tmg_stop(155)
1633
1634 end subroutine transpose2d_LevtoLon
1635
1636 !--------------------------------------------------------------------------
1637 ! lst_Transpose2d_LevToLon_kij_mpitypes
1638 !--------------------------------------------------------------------------
1639 subroutine transpose2d_LevToLon_kij_mpitypes(gd_out,gd_in,nk,lst)
1640 implicit none
1641
1642 ! Arguments:
1643 type(struct_lst), intent(in) :: lst
1644 integer, intent(in) :: nk
1645 real(8), intent(out) :: gd_out(nk,lst%myLonBeg:lst%myLonEnd, lst%myLatBeg:lst%myLatEnd)
1646 real(8), intent(in) :: gd_in(lst%myLevBeg:lst%myLevEnd,lst%ni, lst%myLatBeg:lst%myLatEnd)
1647
1648 ! Locals:
1649 integer :: nsize, ierr
1650
1651 if (verbose) write(*,*) 'Entering transpose2d_LevToLon_kij'
1652 call rpn_comm_barrier("GRID",ierr)
1653
1654 call utl_tmg_start(155,'low-level--lst_transpose_LEVtoLON')
1655
1656 nsize = lst%lonPerPE*lst%maxLevCount*lst%latPerPE
1657 if (mmpi_npex > 1) then
1658 call mpi_alltoall(gd_in, 1, lst%sendType_LevToLon, &
1659 gd_out, 1, lst%recvType_LevToLon, mmpi_comm_EW, ierr)
1660 else
1661 gd_out(:,:,:) = gd_in(:,:,:)
1662 end if
1663
1664 call utl_tmg_stop(155)
1665
1666 end subroutine transpose2d_LevtoLon_kij_mpitypes
1667
1668 !--------------------------------------------------------------------------
1669 ! lst_Transpose2d_LevToLon_kij
1670 !--------------------------------------------------------------------------
1671 subroutine transpose2d_LevToLon_kij(gd_out,gd_in,nk,lst)
1672 implicit none
1673
1674 ! Arguments:
1675 type(struct_lst), intent(in) :: lst
1676 integer, intent(in) :: nk
1677 real(8), intent(out) :: gd_out(nk,lst%myLonBeg:lst%myLonEnd, lst%myLatBeg:lst%myLatEnd)
1678 real(8), intent(in) :: gd_in(lst%myLevBeg:lst%myLevEnd,lst%ni, lst%myLatBeg:lst%myLatEnd)
1679
1680 ! Locals:
1681 real(8) :: gd_send(lst%maxLevCount, lst%lonPerPEmax, lst%latPerPEmax, mmpi_npex)
1682 real(8) :: gd_recv(lst%maxLevCount, lst%lonPerPEmax, lst%latPerPEmax, mmpi_npex)
1683 integer :: yourid, nsize, ierr, levIndex, levIndex2, latIndex, latIndex2, lonIndex, lonIndex2
1684
1685 if (verbose) write(*,*) 'Entering transpose2d_LevToLon_kij'
1686 call rpn_comm_barrier("GRID",ierr)
1687
1688 call utl_tmg_start(155,'low-level--lst_transpose_LEVtoLON')
1689
1690 !$OMP PARALLEL DO PRIVATE(yourid,levIndex,levIndex2,lonIndex,lonIndex2,latIndex,latIndex2)
1691 do yourid = 0, (mmpi_npex-1)
1692 gd_send(:,:,:,yourid+1) = 0.0d0
1693 do latIndex = lst%myLatBeg, lst%myLatEnd
1694 latIndex2 = latIndex - lst%myLatBeg + 1
1695 do lonIndex = lst%allLonBeg(yourid+1), lst%allLonEnd(yourid+1)
1696 lonIndex2 = lonIndex - lst%allLonBeg(yourid+1) + 1
1697 do levIndex = lst%myLevBeg, lst%myLevEnd
1698 levIndex2 = levIndex - lst%myLevBeg + 1
1699 gd_send(levIndex2,lonIndex2,latIndex2,yourid+1) = gd_in(levIndex,lonIndex,latIndex)
1700 end do
1701 end do
1702 end do
1703 end do
1704 !$OMP END PARALLEL DO
1705
1706 nsize = lst%lonPerPEmax * lst%maxLevCount * lst%latPerPEmax
1707 if (mmpi_npex > 1) then
1708 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
1709 gd_recv,nsize,"mpi_double_precision","EW",ierr)
1710 else
1711 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
1712 end if
1713
1714 !$OMP PARALLEL DO PRIVATE(yourid,latIndex,latIndex2,levIndex,levIndex2,lonIndex,lonIndex2)
1715 do yourid = 0, (mmpi_npex-1)
1716 do latIndex = lst%myLatBeg, lst%myLatEnd
1717 latIndex2 = latIndex - lst%myLatBeg + 1
1718 do lonIndex = lst%myLonBeg, lst%myLonEnd
1719 lonIndex2 = lonIndex - lst%myLonBeg + 1
1720 do levIndex=lst%allLevBeg(yourid+1),lst%allLevEnd(yourid+1)
1721 levIndex2=levIndex-lst%allLevBeg(yourid+1)+1
1722 gd_out(levIndex,lonIndex,latIndex) = gd_recv(levIndex2,lonIndex2,latIndex2,yourid+1)
1723 end do
1724 end do
1725 end do
1726 end do
1727 !$OMP END PARALLEL DO
1728
1729 call utl_tmg_stop(155)
1730
1731 end subroutine transpose2d_LevtoLon_kij
1732
1733 !--------------------------------------------------------------------------
1734 ! lst_Transpose2d_LatToM
1735 !--------------------------------------------------------------------------
1736 subroutine transpose2d_LatToM(gd_out, gd_in, lst)
1737 implicit none
1738
1739 ! Arguments:
1740 type(struct_lst), intent(in) :: lst
1741 real(8), intent(out) :: gd_out(2*lst%mymCount,lst%nj+lst%njp,lst%myLevBeg:lst%myLevEnd)
1742 real(8), intent(in) :: gd_in (lst%ni+lst%nip,lst%latPerPE,lst%myLevBeg:lst%myLevEnd)
1743
1744 ! Locals:
1745 real(8) :: gd_recv(lst%maxmActiveCount,2,lst%latPerPEmax, lst%maxLevCount, mmpi_npey)
1746 real(8) :: gd_send(lst%maxmActiveCount,2,lst%latPerPEmax, lst%maxLevCount, mmpi_npey)
1747 integer :: yourid, mIndex, icount, nsize, ierr, levIndex, levIndex2, latIndex, latIndex2
1748
1749 if (verbose) write(*,*) 'Entering transpose2d_LatToM'
1750 call rpn_comm_barrier("GRID",ierr)
1751
1752 call utl_tmg_start(154,'low-level--lst_transpose_MtoLAT')
1753
1754 !$OMP PARALLEL DO PRIVATE(yourid,latIndex,levIndex,levIndex2,icount,mIndex)
1755 do yourid = 0, (mmpi_npey-1)
1756 do levIndex = lst%myLevBeg, lst%myLevEnd
1757 levIndex2 = levIndex - lst%myLevBeg + 1
1758 gd_send(:,:,:,levIndex2,yourid+1) = 0.d0
1759 do latIndex = 1, lst%latPerPE
1760 icount = 0
1761 do mIndex = lst%allmBeg(yourid+1), lst%allmEnd(yourid+1), lst%allmSkip(yourid+1)
1762 if (lst%KfromMNglb(mIndex,0) /= -1) then
1763 icount = icount + 1
1764 gd_send(icount,1,latIndex,levIndex2,yourid+1) = gd_in(2*mIndex+1, latIndex, levIndex)
1765 gd_send(icount,2,latIndex,levIndex2,yourid+1) = gd_in(2*mIndex+2, latIndex, levIndex)
1766 end if
1767 end do
1768 end do
1769 end do
1770 end do
1771 !$OMP END PARALLEL DO
1772
1773 nsize = lst%maxmActiveCount * 2 * lst%maxLevCount * lst%latPerPEmax
1774 if (mmpi_npey > 1) then
1775 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
1776 gd_recv,nsize,"mpi_double_precision","NS",ierr)
1777 else
1778 gd_recv(:,:,:,:,1) = gd_send(:,:,:,:,1)
1779 end if
1780
1781 !$OMP PARALLEL DO PRIVATE(yourid,latIndex,latIndex2,levIndex,levIndex2,icount,mIndex)
1782 do yourid = 0, (mmpi_npey-1)
1783 do levIndex = lst%myLevBeg, lst%myLevEnd
1784 levIndex2 = levIndex - lst%myLevBeg + 1
1785 do latIndex = lst%allLatBeg(yourid+1), lst%allLatEnd(yourid+1)
1786 latIndex2 = latIndex - lst%allLatBeg(yourid+1) + 1
1787 icount = 0
1788 do mIndex = lst%mymBeg, lst%mymEnd, lst%mymSkip
1789 if (lst%KfromMNglb(mIndex,0) /= -1) then
1790 icount = icount + 1
1791 gd_out(2*lst%mymIndex(mIndex)-1,latIndex,levIndex) = gd_recv(icount,1,latIndex2,levIndex2,yourid+1)
1792 gd_out(2*lst%mymIndex(mIndex) ,latIndex,levIndex) = gd_recv(icount,2,latIndex2,levIndex2,yourid+1)
1793 else
1794 gd_out(2*lst%mymIndex(mIndex)-1,latIndex,levIndex) = 0.d0
1795 gd_out(2*lst%mymIndex(mIndex) ,latIndex,levIndex) = 0.d0
1796 end if
1797 end do
1798 end do
1799 end do
1800 end do
1801 !$OMP END PARALLEL DO
1802
1803 call utl_tmg_stop(154)
1804
1805 end subroutine transpose2d_LatToM
1806
1807 !--------------------------------------------------------------------------
1808 ! lst_Transpose2d_LatToM_kij
1809 !--------------------------------------------------------------------------
1810 subroutine transpose2d_LatToM_kij(gd_out, gd_in, lst)
1811 implicit none
1812
1813 ! Arguments:
1814 type(struct_lst), intent(in) :: lst
1815 real(8), intent(out) :: gd_out(lst%myLevBeg:lst%myLevEnd,2*lst%mymCount,lst%nj+lst%njp )
1816 real(8), intent(in) :: gd_in (lst%myLevBeg:lst%myLevEnd,lst%ni+lst%nip,lst%latPerPE)
1817
1818 ! Locals:
1819 real(8) :: gd_recv(lst%maxLevCount,lst%maxmActiveCount,2,lst%latPerPEmax, mmpi_npey)
1820 real(8) :: gd_send(lst%maxLevCount,lst%maxmActiveCount,2,lst%latPerPEmax, mmpi_npey)
1821 integer :: yourid, mIndex, icount, nsize, ierr, levIndex, levIndex2, latIndex, latIndex2
1822
1823 if (verbose) write(*,*) 'Entering transpose2d_LatToM_kij'
1824 call rpn_comm_barrier("GRID",ierr)
1825
1826 call utl_tmg_start(154,'low-level--lst_transpose_MtoLAT')
1827
1828 !$OMP PARALLEL DO PRIVATE(yourid,latIndex,levIndex,levIndex2,icount,mIndex)
1829 do yourid = 0, (mmpi_npey-1)
1830 gd_send(:,:,:,:,yourid+1) = 0.d0
1831 do latIndex = 1, lst%latPerPE
1832 icount = 0
1833 do mIndex = lst%allmBeg(yourid+1), lst%allmEnd(yourid+1), lst%allmSkip(yourid+1)
1834 if (lst%KfromMNglb(mIndex,0) /= -1) then
1835 icount = icount + 1
1836 do levIndex = lst%myLevBeg, lst%myLevEnd
1837 levIndex2 = levIndex - lst%myLevBeg + 1
1838 gd_send(levIndex2,icount,1,latIndex,yourid+1) = gd_in(levIndex,2*mIndex+1, latIndex)
1839 gd_send(levIndex2,icount,2,latIndex,yourid+1) = gd_in(levIndex,2*mIndex+2, latIndex)
1840 end do
1841 end if
1842 end do
1843 end do
1844 end do
1845 !$OMP END PARALLEL DO
1846
1847 nsize = lst%maxmActiveCount * 2 * lst%maxLevCount * lst%latPerPEmax
1848 if (mmpi_npey > 1) then
1849 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
1850 gd_recv,nsize,"mpi_double_precision","NS",ierr)
1851 else
1852 gd_recv(:,:,:,:,1) = gd_send(:,:,:,:,1)
1853 end if
1854
1855 !$OMP PARALLEL DO PRIVATE(yourid,latIndex,latIndex2,levIndex,levIndex2,icount,mIndex)
1856 do yourid = 0, (mmpi_npey-1)
1857 do latIndex = lst%allLatBeg(yourid+1), lst%allLatEnd(yourid+1)
1858 latIndex2 = latIndex - lst%allLatBeg(yourid+1) + 1
1859 icount = 0
1860 do mIndex = lst%mymBeg, lst%mymEnd, lst%mymSkip
1861 if (lst%KfromMNglb(mIndex,0) /= -1) then
1862 icount = icount + 1
1863 do levIndex = lst%myLevBeg, lst%myLevEnd
1864 levIndex2 = levIndex - lst%myLevBeg + 1
1865 gd_out(levIndex,2*lst%mymIndex(mIndex)-1,latIndex) = gd_recv(levIndex2,icount,1,latIndex2,yourid+1)
1866 gd_out(levIndex,2*lst%mymIndex(mIndex) ,latIndex) = gd_recv(levIndex2,icount,2,latIndex2,yourid+1)
1867 end do
1868 else
1869 gd_out(:,2*lst%mymIndex(mIndex)-1,latIndex) = 0.d0
1870 gd_out(:,2*lst%mymIndex(mIndex) ,latIndex) = 0.d0
1871 end if
1872 end do
1873 end do
1874 end do
1875 !$OMP END PARALLEL DO
1876
1877 call utl_tmg_stop(154)
1878
1879 end subroutine transpose2d_LatToM_kij
1880
1881 !--------------------------------------------------------------------------
1882 ! Transpose2d_MToLat
1883 !--------------------------------------------------------------------------
1884 subroutine transpose2d_MtoLat(gd_out, gd_in, lst)
1885 implicit none
1886
1887 ! Arguments:
1888 type(struct_lst), intent(in) :: lst
1889 real(8), intent(in) :: gd_in (2*lst%mymCount,lst%nj+lst%njp,lst%myLevBeg:lst%myLevEnd)
1890 real(8), intent(out) :: gd_out(lst%ni+lst%nip,lst%latPerPE,lst%myLevBeg:lst%myLevEnd)
1891
1892 ! Locals:
1893 real(8) :: gd_recv(lst%maxmActiveCount,2,lst%latPerPEmax,lst%maxLevCount, mmpi_npey)
1894 real(8) :: gd_send(lst%maxmActiveCount,2,lst%latPerPEmax,lst%maxLevCount, mmpi_npey)
1895 integer :: yourid, mIndex, icount, nsize, ierr, levIndex, levIndex2, latIndex, latIndex2
1896
1897 if (verbose) write(*,*) 'Entering transpose2d_MToLat'
1898 call rpn_comm_barrier("GRID",ierr)
1899
1900 call utl_tmg_start(154,'low-level--lst_transpose_MtoLAT')
1901
1902 !$OMP PARALLEL DO PRIVATE(yourid,latIndex,latIndex2,levIndex,levIndex2,icount,mIndex)
1903 do yourid = 0, (mmpi_npey-1)
1904 do levIndex = lst%myLevBeg, lst%myLevEnd
1905 levIndex2 = levIndex - lst%myLevBeg + 1
1906 gd_send(:,:,:,levIndex2,yourid+1) = 0.d0
1907 do latIndex = lst%allLatBeg(yourid+1), lst%allLatEnd(yourid+1)
1908 latIndex2 = latIndex - lst%allLatBeg(yourid+1) + 1
1909 icount = 0
1910 do mIndex = lst%mymBeg, lst%mymEnd, lst%mymSkip
1911 if (lst%KfromMNglb(mIndex,0) /= -1) then
1912 icount = icount+1
1913 gd_send(icount,1,latIndex2,levIndex2,yourid+1) = gd_in(2*lst%mymIndex(mIndex)-1,latIndex,levIndex)
1914 gd_send(icount,2,latIndex2,levIndex2,yourid+1) = gd_in(2*lst%mymIndex(mIndex) ,latIndex,levIndex)
1915 end if
1916 end do
1917 end do
1918 end do
1919 end do
1920 !$OMP END PARALLEL DO
1921
1922 nsize = lst%maxmActiveCount * 2 * lst%maxLevCount * lst%latPerPEmax
1923 if (mmpi_npey > 1) then
1924 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
1925 gd_recv,nsize,"mpi_double_precision","NS",ierr)
1926 else
1927 gd_recv(:,:,:,:,1) = gd_send(:,:,:,:,1)
1928 end if
1929
1930 !$OMP PARALLEL DO PRIVATE(yourid,latIndex,levIndex,levIndex2,icount,mIndex)
1931 do yourid = 0, (mmpi_npey-1)
1932 do levIndex = lst%myLevBeg, lst%myLevEnd
1933 levIndex2 = levIndex - lst%myLevBeg + 1
1934 do latIndex = 1, lst%latPerPE
1935 icount = 0
1936 do mIndex = lst%allmBeg(yourid+1), lst%allmEnd(yourid+1), lst%allmSkip(yourid+1)
1937 if (lst%KfromMNglb(mIndex,0) /= -1) then
1938 icount = icount+1
1939 gd_out(2*mIndex+1,latIndex,levIndex) = gd_recv(icount,1,latIndex,levIndex2,yourid+1)
1940 gd_out(2*mIndex+2,latIndex,levIndex) = gd_recv(icount,2,latIndex,levIndex2,yourid+1)
1941 else
1942 gd_out(2*mIndex+1,latIndex,levIndex) = 0.d0
1943 gd_out(2*mIndex+2,latIndex,levIndex) = 0.d0
1944 end if
1945 end do
1946 end do
1947 end do
1948 end do
1949 !$OMP END PARALLEL DO
1950
1951 call utl_tmg_stop(154)
1952
1953 end subroutine transpose2d_MtoLat
1954
1955 !--------------------------------------------------------------------------
1956 ! Transpose2d_MToLat_kij
1957 !--------------------------------------------------------------------------
1958 subroutine transpose2d_MtoLat_kij(gd_out, gd_in, lst)
1959 implicit none
1960
1961 ! Arguments:
1962 type(struct_lst), intent(in) :: lst
1963 real(8), intent(in) :: gd_in (lst%myLevBeg:lst%myLevEnd,2*lst%mymCount,lst%nj+lst%njp)
1964 real(8), intent(out) :: gd_out(lst%myLevBeg:lst%myLevEnd,lst%ni+lst%nip,lst%latPerPE)
1965
1966 ! Locals:
1967 real(8) :: gd_recv(lst%maxLevCount,lst%maxmActiveCount,2,lst%latPerPEmax, mmpi_npey)
1968 real(8) :: gd_send(lst%maxLevCount,lst%maxmActiveCount,2,lst%latPerPEmax, mmpi_npey)
1969 integer :: yourid, mIndex, icount, nsize, ierr, levIndex, levIndex2, latIndex, latIndex2
1970
1971 if (verbose) write(*,*) 'Entering transpose2d_MToLat_kij'
1972 call rpn_comm_barrier("GRID",ierr)
1973
1974 call utl_tmg_start(154,'low-level--lst_transpose_MtoLAT')
1975
1976 !$OMP PARALLEL DO PRIVATE(yourid,latIndex,latIndex2,levIndex,levIndex2,icount,mIndex)
1977 do yourid = 0, (mmpi_npey-1)
1978 gd_send(:,:,:,:,yourid+1) = 0.d0
1979 do latIndex = lst%allLatBeg(yourid+1), lst%allLatEnd(yourid+1)
1980 latIndex2 = latIndex - lst%allLatBeg(yourid+1) + 1
1981 icount = 0
1982 do mIndex = lst%mymBeg, lst%mymEnd, lst%mymSkip
1983 if (lst%KfromMNglb(mIndex,0) /= -1) then
1984 icount = icount+1
1985 do levIndex = lst%myLevBeg, lst%myLevEnd
1986 levIndex2 = levIndex - lst%myLevBeg + 1
1987 gd_send(levIndex2,icount,1,latIndex2,yourid+1) = gd_in(levIndex,2*lst%mymIndex(mIndex)-1,latIndex)
1988 gd_send(levIndex2,icount,2,latIndex2,yourid+1) = gd_in(levIndex,2*lst%mymIndex(mIndex) ,latIndex)
1989 end do
1990 end if
1991 end do
1992 end do
1993 end do
1994 !$OMP END PARALLEL DO
1995
1996 nsize = lst%maxmActiveCount * 2 * lst%maxLevCount * lst%latPerPEmax
1997 if (mmpi_npey > 1) then
1998 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
1999 gd_recv,nsize,"mpi_double_precision","NS",ierr)
2000 else
2001 gd_recv(:,:,:,:,1) = gd_send(:,:,:,:,1)
2002 end if
2003
2004 !$OMP PARALLEL DO PRIVATE(yourid,latIndex,levIndex,levIndex2,icount,mIndex)
2005 do yourid = 0, (mmpi_npey-1)
2006 do latIndex = 1, lst%latPerPE
2007 icount = 0
2008 do mIndex = lst%allmBeg(yourid+1), lst%allmEnd(yourid+1), lst%allmSkip(yourid+1)
2009 if (lst%KfromMNglb(mIndex,0) /= -1) then
2010 icount = icount+1
2011 do levIndex = lst%myLevBeg, lst%myLevEnd
2012 levIndex2 = levIndex - lst%myLevBeg + 1
2013 gd_out(levIndex,2*mIndex+1,latIndex) = gd_recv(levIndex2,icount,1,latIndex,yourid+1)
2014 gd_out(levIndex,2*mIndex+2,latIndex) = gd_recv(levIndex2,icount,2,latIndex,yourid+1)
2015 end do
2016 else
2017 gd_out(:,2*mIndex+1,latIndex) = 0.d0
2018 gd_out(:,2*mIndex+2,latIndex) = 0.d0
2019 end if
2020 end do
2021 end do
2022 end do
2023 !$OMP END PARALLEL DO
2024
2025 call utl_tmg_stop(154)
2026
2027 end subroutine transpose2d_MtoLat_kij
2028
2029 !--------------------------------------------------------------------------
2030 ! lst_Transpose2d_LevToN
2031 !--------------------------------------------------------------------------
2032 subroutine transpose2d_LevToN(SpectralStateVar, gd_in, nk, lst)
2033 implicit none
2034
2035 ! Arguments:
2036 type(struct_lst), intent(in) :: lst
2037 integer, intent(in) :: nk
2038 real(8), intent(out) :: SpectralStateVar(lst%nla,lst%nphase,nk)
2039 real(8), intent(in) :: gd_in (2*lst%mymCount, lst%nj+lst%njp, lst%myLevBeg:lst%myLevEnd)
2040
2041 ! Locals:
2042 real(8) :: gd_send(lst%maxnla, 4, lst%maxLevCount, mmpi_npex)
2043 real(8) :: gd_recv(lst%maxnla, 4, lst%maxLevCount, mmpi_npex)
2044 integer :: yourid, nsize, ierr, levIndex, levIndex2, nIndex, mIndex, icount
2045
2046 if (verbose) write(*,*) 'Entering transpose2d_LevToN'
2047 call rpn_comm_barrier("GRID",ierr)
2048
2049 call utl_tmg_start(153,'low-level--lst_transpose_NtoLEV')
2050
2051 !$OMP PARALLEL DO PRIVATE(yourid,mIndex,levIndex,levIndex2,nIndex,icount)
2052 do yourid = 0, (mmpi_npex-1)
2053 do levIndex = lst%myLevBeg, lst%myLevEnd
2054 levIndex2 = levIndex - lst%myLevBeg + 1
2055 gd_send(:,:,levIndex2,yourid+1) = 0.d0
2056 icount = 0
2057 do nIndex = lst%allnBeg(yourid+1), lst%allnEnd(yourid+1), lst%allnSkip(yourid+1)
2058 do mIndex = lst%mymBeg, lst%mymEnd, lst%mymSkip
2059 if (lst%KfromMNglb(mIndex,nIndex) /= -1) then
2060 icount = icount + 1
2061 gd_send(icount,1,levIndex2,yourid+1) = gd_in(2*lst%mymIndex(mIndex)-1,2*nIndex+1,levIndex)
2062 gd_send(icount,2,levIndex2,yourid+1) = gd_in(2*lst%mymIndex(mIndex)-1,2*nIndex+2,levIndex)
2063 gd_send(icount,3,levIndex2,yourid+1) = gd_in(2*lst%mymIndex(mIndex) ,2*nIndex+1,levIndex)
2064 gd_send(icount,4,levIndex2,yourid+1) = gd_in(2*lst%mymIndex(mIndex) ,2*nIndex+2,levIndex)
2065 end if
2066 end do
2067 end do
2068 end do
2069 end do
2070 !$OMP END PARALLEL DO
2071
2072 nsize = lst%maxnla * 4 * lst%maxLevCount
2073 if (mmpi_npex > 1) then
2074 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
2075 gd_recv,nsize,"mpi_double_precision","EW",ierr)
2076 else
2077 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
2078 end if
2079
2080 !$OMP PARALLEL DO PRIVATE(yourid,levIndex,levIndex2)
2081 do yourid = 0, (mmpi_npex-1)
2082 do levIndex = lst%allLevBeg(yourid+1), lst%allLevEnd(yourid+1)
2083 levIndex2 = levIndex - lst%allLevBeg(yourid+1) + 1
2084 SpectralStateVar(:,:,levIndex) = gd_recv(1:lst%nla,:,levIndex2,yourid+1)
2085 end do
2086 end do
2087 !$OMP END PARALLEL DO
2088
2089 call utl_tmg_stop(153)
2090
2091 end subroutine transpose2d_LevToN
2092
2093 !--------------------------------------------------------------------------
2094 ! lst_Transpose2d_LevToN_kij
2095 !--------------------------------------------------------------------------
2096 subroutine transpose2d_LevToN_kij(SpectralStateVar, gd_in, nk, lst)
2097 implicit none
2098
2099 ! Arguments:
2100 type(struct_lst), intent(in) :: lst
2101 integer, intent(in) :: nk
2102 real(8), intent(out) :: SpectralStateVar(lst%nla,lst%nphase,nk)
2103 real(8), intent(in) :: gd_in (lst%myLevBeg:lst%myLevEnd,2*lst%mymCount, lst%nj+lst%njp)
2104
2105 ! Locals:
2106 real(8) :: gd_send(lst%maxnla, 4, lst%maxLevCount, mmpi_npex)
2107 real(8) :: gd_recv(lst%maxnla, 4, lst%maxLevCount, mmpi_npex)
2108 integer :: yourid, nsize, ierr, levIndex, levIndex2, nIndex, mIndex, icount
2109
2110 if (verbose) write(*,*) 'Entering transpose2d_LevToN_kij'
2111 call rpn_comm_barrier("GRID",ierr)
2112
2113 call utl_tmg_start(153,'low-level--lst_transpose_NtoLEV')
2114
2115 !$OMP PARALLEL DO PRIVATE(yourid,mIndex,levIndex,levIndex2,nIndex,icount)
2116 do yourid = 0, (mmpi_npex-1)
2117 do levIndex = lst%myLevBeg, lst%myLevEnd
2118 levIndex2 = levIndex - lst%myLevBeg + 1
2119 gd_send(:,:,levIndex2,yourid+1) = 0.d0
2120 icount = 0
2121 do nIndex = lst%allnBeg(yourid+1), lst%allnEnd(yourid+1), lst%allnSkip(yourid+1)
2122 do mIndex = lst%mymBeg, lst%mymEnd, lst%mymSkip
2123 if (lst%KfromMNglb(mIndex,nIndex) /= -1) then
2124 icount = icount + 1
2125 gd_send(icount,1,levIndex2,yourid+1) = gd_in(levIndex,2*lst%mymIndex(mIndex)-1,2*nIndex+1)
2126 gd_send(icount,2,levIndex2,yourid+1) = gd_in(levIndex,2*lst%mymIndex(mIndex)-1,2*nIndex+2)
2127 gd_send(icount,3,levIndex2,yourid+1) = gd_in(levIndex,2*lst%mymIndex(mIndex) ,2*nIndex+1)
2128 gd_send(icount,4,levIndex2,yourid+1) = gd_in(levIndex,2*lst%mymIndex(mIndex) ,2*nIndex+2)
2129 end if
2130 end do
2131 end do
2132 end do
2133 end do
2134 !$OMP END PARALLEL DO
2135
2136 nsize = lst%maxnla * 4 * lst%maxLevCount
2137 if (mmpi_npex > 1) then
2138 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
2139 gd_recv,nsize,"mpi_double_precision","EW",ierr)
2140 else
2141 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
2142 end if
2143
2144 !$OMP PARALLEL DO PRIVATE(yourid,levIndex,levIndex2)
2145 do yourid = 0, (mmpi_npex-1)
2146 do levIndex = lst%allLevBeg(yourid+1), lst%allLevEnd(yourid+1)
2147 levIndex2 = levIndex - lst%allLevBeg(yourid+1) + 1
2148 SpectralStateVar(:,:,levIndex) = gd_recv(1:lst%nla,:,levIndex2,yourid+1)
2149 end do
2150 end do
2151 !$OMP END PARALLEL DO
2152
2153 call utl_tmg_stop(153)
2154
2155 end subroutine transpose2d_LevToN_kij
2156
2157 !--------------------------------------------------------------------------
2158 ! lst_Transpose2d_NToLev
2159 !--------------------------------------------------------------------------
2160 subroutine transpose2d_NToLev(gd_out, SpectralStateVar, nk, lst)
2161 implicit none
2162
2163 ! Arguments:
2164 type(struct_lst), intent(in) :: lst
2165 integer, intent(in) :: nk
2166 real(8), intent(in) :: SpectralStateVar(lst%nla,lst%nphase,nk)
2167 real(8), intent(out) :: gd_out(2*lst%mymCount, lst%nj+lst%njp, lst%myLevBeg:lst%myLevEnd)
2168
2169 ! Locals:
2170 real(8) :: gd_send(lst%maxnla, 4, lst%maxLevCount, mmpi_npex)
2171 real(8) :: gd_recv(lst%maxnla, 4, lst%maxLevCount, mmpi_npex)
2172 integer :: yourid, nsize, ierr, levIndex, levIndex2, nIndex, mIndex, icount
2173
2174 if (verbose) write(*,*) 'Entering transpose2d_NToLev'
2175 call rpn_comm_barrier("GRID",ierr)
2176
2177 call utl_tmg_start(153,'low-level--lst_transpose_NtoLEV')
2178
2179 !$OMP PARALLEL DO PRIVATE(yourid,levIndex,levIndex2)
2180 do yourid = 0, (mmpi_npex-1)
2181 do levIndex = lst%allLevBeg(yourid+1), lst%allLevEnd(yourid+1)
2182 levIndex2 = levIndex - lst%allLevBeg(yourid+1) + 1
2183 gd_send(1:lst%nla,:,levIndex2,yourid+1) = SpectralStateVar(:,:,levIndex)
2184 if (lst%nla < lst%maxnla) gd_send(lst%nla+1:lst%maxnla,:,levIndex2,yourid+1) = 0.d0
2185 end do
2186 end do
2187 !$OMP END PARALLEL DO
2188
2189 nsize = lst%maxnla * 4* lst%maxLevCount
2190 if (mmpi_npex > 1) then
2191 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
2192 gd_recv,nsize,"mpi_double_precision","EW",ierr)
2193 else
2194 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
2195 end if
2196
2197 !$OMP PARALLEL DO PRIVATE(yourid,nIndex,levIndex,levIndex2,mIndex,icount)
2198 do yourid = 0, (mmpi_npex-1)
2199 do levIndex = lst%myLevBeg, lst%myLevEnd
2200 levIndex2 = levIndex - lst%myLevBeg + 1
2201 icount = 0
2202 do nIndex = lst%allnBeg(yourid+1), lst%allnEnd(yourid+1), lst%allnSkip(yourid+1)
2203 do mIndex = lst%mymBeg, lst%mymEnd, lst%mymSkip
2204 if (lst%KfromMNglb(mIndex,nIndex) /= -1) then
2205 icount = icount + 1
2206 gd_out(2*lst%mymIndex(mIndex)-1,2*nIndex+1,levIndex) = gd_recv(icount,1,levIndex2,yourid+1)
2207 gd_out(2*lst%mymIndex(mIndex)-1,2*nIndex+2,levIndex) = gd_recv(icount,2,levIndex2,yourid+1)
2208 gd_out(2*lst%mymIndex(mIndex) ,2*nIndex+1,levIndex) = gd_recv(icount,3,levIndex2,yourid+1)
2209 gd_out(2*lst%mymIndex(mIndex) ,2*nIndex+2,levIndex) = gd_recv(icount,4,levIndex2,yourid+1)
2210 else
2211 gd_out(2*lst%mymIndex(mIndex)-1,2*nIndex+1,levIndex) = 0.d0
2212 gd_out(2*lst%mymIndex(mIndex)-1,2*nIndex+2,levIndex) = 0.d0
2213 gd_out(2*lst%mymIndex(mIndex) ,2*nIndex+1,levIndex) = 0.d0
2214 gd_out(2*lst%mymIndex(mIndex) ,2*nIndex+2,levIndex) = 0.d0
2215 end if
2216 end do
2217 end do
2218 end do
2219 end do
2220 !$OMP END PARALLEL DO
2221
2222 call utl_tmg_stop(153)
2223
2224 end subroutine transpose2d_NToLev
2225
2226 !--------------------------------------------------------------------------
2227 ! lst_Transpose2d_NToLev_kij
2228 !--------------------------------------------------------------------------
2229 subroutine transpose2d_NToLev_kij(gd_out, SpectralStateVar, nk, lst)
2230 implicit none
2231
2232 ! Arguments:
2233 type(struct_lst), intent(in) :: lst
2234 integer, intent(in) :: nk
2235 real(8), intent(in) :: SpectralStateVar(lst%nla,lst%nphase,nk)
2236 real(8), intent(out) :: gd_out(lst%myLevBeg:lst%myLevEnd,2*lst%mymCount,lst%nj+lst%njp)
2237
2238 ! Locals:
2239 real(8) :: gd_send(lst%maxnla, 4, lst%maxLevCount, mmpi_npex)
2240 real(8) :: gd_recv(lst%maxnla, 4, lst%maxLevCount, mmpi_npex)
2241 integer :: yourid, nsize, ierr, levIndex, levIndex2, nIndex, mIndex, icount
2242
2243 if (verbose) write(*,*) 'Entering transpose2d_NToLev_kij'
2244 call rpn_comm_barrier("GRID",ierr)
2245
2246 call utl_tmg_start(153,'low-level--lst_transpose_NtoLEV')
2247
2248 !$OMP PARALLEL DO PRIVATE(yourid,levIndex,levIndex2)
2249 do yourid = 0, (mmpi_npex-1)
2250 do levIndex = lst%allLevBeg(yourid+1), lst%allLevEnd(yourid+1)
2251 levIndex2 = levIndex - lst%allLevBeg(yourid+1) + 1
2252 gd_send(1:lst%nla,:,levIndex2,yourid+1) = SpectralStateVar(:,:,levIndex)
2253 if (lst%nla < lst%maxnla) gd_send(lst%nla+1:lst%maxnla,:,levIndex2,yourid+1) = 0.d0
2254 end do
2255 end do
2256 !$OMP END PARALLEL DO
2257
2258 nsize = lst%maxnla * 4 * lst%maxLevCount
2259 if (mmpi_npex > 1) then
2260 call rpn_comm_alltoall(gd_send,nsize,"mpi_double_precision", &
2261 gd_recv,nsize,"mpi_double_precision","EW",ierr)
2262 else
2263 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
2264 end if
2265
2266 !$OMP PARALLEL DO PRIVATE(yourid,nIndex,levIndex,levIndex2,mIndex,icount)
2267 do yourid = 0, (mmpi_npex-1)
2268 do levIndex = lst%myLevBeg, lst%myLevEnd
2269 levIndex2 = levIndex - lst%myLevBeg + 1
2270 icount = 0
2271 do nIndex = lst%allnBeg(yourid+1), lst%allnEnd(yourid+1), lst%allnSkip(yourid+1)
2272 do mIndex = lst%mymBeg, lst%mymEnd, lst%mymSkip
2273 if (lst%KfromMNglb(mIndex,nIndex) /= -1) then
2274 icount = icount + 1
2275 gd_out(levIndex,2*lst%mymIndex(mIndex)-1,2*nIndex+1) = gd_recv(icount,1,levIndex2,yourid+1)
2276 gd_out(levIndex,2*lst%mymIndex(mIndex)-1,2*nIndex+2) = gd_recv(icount,2,levIndex2,yourid+1)
2277 gd_out(levIndex,2*lst%mymIndex(mIndex) ,2*nIndex+1) = gd_recv(icount,3,levIndex2,yourid+1)
2278 gd_out(levIndex,2*lst%mymIndex(mIndex) ,2*nIndex+2) = gd_recv(icount,4,levIndex2,yourid+1)
2279 else
2280 gd_out(levIndex,2*lst%mymIndex(mIndex)-1,2*nIndex+1) = 0.d0
2281 gd_out(levIndex,2*lst%mymIndex(mIndex)-1,2*nIndex+2) = 0.d0
2282 gd_out(levIndex,2*lst%mymIndex(mIndex) ,2*nIndex+1) = 0.d0
2283 gd_out(levIndex,2*lst%mymIndex(mIndex) ,2*nIndex+2) = 0.d0
2284 end if
2285 end do
2286 end do
2287 end do
2288 end do
2289 !$OMP END PARALLEL DO
2290
2291 call utl_tmg_stop(153)
2292
2293 end subroutine transpose2d_NToLev_kij
2294
2295 !--------------------------------------------------------------------------
2296 ! lst_ReshapeTrunc
2297 !--------------------------------------------------------------------------
2298 subroutine lst_ReshapeTrunc(SpectralStateRpn, SpectralStateVar, &
2299 Direction, kStart, kEnd, lst)
2300 implicit none
2301
2302 ! Arguments:
2303 type(struct_lst), intent(in) :: lst
2304 integer, intent(in) :: kStart
2305 integer, intent(in) :: kEnd
2306 character(len=*), intent(in) :: Direction ! ToVAR or ToRPN
2307 real(8), intent(inout) :: SpectralStateRpn(2*lst%mymCount,2*lst%mynCount,kStart:kEnd)
2308 real(8), intent(inout) :: SpectralStateVar(lst%nla ,lst%nphase,kStart:kEnd)
2309
2310 ! Locals:
2311 integer k, m, n, ila
2312
2313 if (verbose) write(*,*) 'Entering lst_ReshapeTrunc'
2314
2315 select case (trim(Direction))
2316 case ('ToVAR')
2317 ! Truncation (if applicable) will be applied here
2318 !$OMP PARALLEL DO PRIVATE (n,m,ila,k)
2319 do n = lst%mynBeg, lst%mynEnd, lst%mynSkip
2320 do m = lst%mymBeg, lst%mymEnd, lst%mymSkip
2321 ila = lst%nla_Index(m,n)
2322 if (ila /= -1) then
2323 do k = kStart, kEnd
2324 SpectralStateVar(ila,1,k) = SpectralStateRpn(2*lst%mymIndex(m)-1,2*lst%mynIndex(n)-1,k)
2325 SpectralStateVar(ila,2,k) = SpectralStateRpn(2*lst%mymIndex(m)-1,2*lst%mynIndex(n) ,k)
2326 SpectralStateVar(ila,3,k) = SpectralStateRpn(2*lst%mymIndex(m) ,2*lst%mynIndex(n)-1,k)
2327 SpectralStateVar(ila,4,k) = SpectralStateRpn(2*lst%mymIndex(m) ,2*lst%mynIndex(n), k)
2328 end do
2329 end if
2330 end do
2331 end do
2332 !$OMP END PARALLEL DO
2333
2334 case ('ToRPN')
2335 SpectralStateRpn(:,:,:) = 0.0d0
2336 !$OMP PARALLEL DO PRIVATE (n,m,ila,k)
2337 do n = lst%mynBeg, lst%mynEnd, lst%mynSkip
2338 do m = lst%mymBeg, lst%mymEnd, lst%mymSkip
2339 ila = lst%nla_Index(m,n)
2340 if (ila /= -1) then
2341 do k = kStart, kEnd
2342 SpectralStateRpn(2*lst%mymIndex(m)-1,2*lst%mynIndex(n)-1,k) = SpectralStateVar(ila,1,k)
2343 SpectralStateRpn(2*lst%mymIndex(m)-1,2*lst%mynIndex(n) ,k) = SpectralStateVar(ila,2,k)
2344 SpectralStateRpn(2*lst%mymIndex(m) ,2*lst%mynIndex(n)-1,k) = SpectralStateVar(ila,3,k)
2345 SpectralStateRpn(2*lst%mymIndex(m) ,2*lst%mynIndex(n) ,k) = SpectralStateVar(ila,4,k)
2346 end do
2347 end if
2348 end do
2349 end do
2350 !$OMP END PARALLEL DO
2351
2352 case default
2353 write(*,*)
2354 write(*,*) 'lst_ReshapeTrunc: Unknown Direction', trim(Direction)
2355 call utl_abort('lst_ReshapeTrunc')
2356 end select
2357
2358 end subroutine lst_ReshapeTrunc
2359
2360 !--------------------------------------------------------------------------
2361 ! lst_ReshapeTrunc_kij
2362 !--------------------------------------------------------------------------
2363 subroutine lst_ReshapeTrunc_kij(SpectralStateRpn, SpectralStateVar, &
2364 Direction, kStart, kEnd, lst)
2365 implicit none
2366
2367 ! Arguments:
2368 type(struct_lst), intent(in) :: lst
2369 integer, intent(in) :: kStart
2370 integer, intent(in) :: kEnd
2371 character(len=*), intent(in) :: Direction ! ToVAR or ToRPN
2372 real(8), intent(inout) :: SpectralStateRpn(kStart:kEnd,2*lst%mymCount,2*lst%mynCount)
2373 real(8), intent(inout) :: SpectralStateVar(lst%nla ,lst%nphase,kStart:kEnd)
2374
2375 ! Locals:
2376 integer k, m, n, ila
2377
2378 if (verbose) write(*,*) 'Entering lst_ReshapeTrunc_kij'
2379
2380 select case (trim(Direction))
2381 case ('ToVAR')
2382 ! Truncation (if applicable) will be applied here
2383 !$OMP PARALLEL DO PRIVATE (n,m,ila,k)
2384 do n = lst%mynBeg, lst%mynEnd, lst%mynSkip
2385 do m = lst%mymBeg, lst%mymEnd, lst%mymSkip
2386 ila = lst%nla_Index(m,n)
2387 if (ila /= -1) then
2388 do k = kStart, kEnd
2389 SpectralStateVar(ila,1,k) = SpectralStateRpn(k,2*lst%mymIndex(m)-1,2*lst%mynIndex(n)-1)
2390 SpectralStateVar(ila,2,k) = SpectralStateRpn(k,2*lst%mymIndex(m)-1,2*lst%mynIndex(n) )
2391 SpectralStateVar(ila,3,k) = SpectralStateRpn(k,2*lst%mymIndex(m) ,2*lst%mynIndex(n)-1)
2392 SpectralStateVar(ila,4,k) = SpectralStateRpn(k,2*lst%mymIndex(m) ,2*lst%mynIndex(n) )
2393 end do
2394 end if
2395 end do
2396 end do
2397 !$OMP END PARALLEL DO
2398
2399 case ('ToRPN')
2400 SpectralStateRpn(:,:,:) = 0.0d0
2401 !$OMP PARALLEL DO PRIVATE (n,m,ila,k)
2402 do n = lst%mynBeg, lst%mynEnd, lst%mynSkip
2403 do m = lst%mymBeg, lst%mymEnd, lst%mymSkip
2404 ila = lst%nla_Index(m,n)
2405 if (ila /= -1) then
2406 do k = kStart, kEnd
2407 SpectralStateRpn(k,2*lst%mymIndex(m)-1,2*lst%mynIndex(n)-1) = SpectralStateVar(ila,1,k)
2408 SpectralStateRpn(k,2*lst%mymIndex(m)-1,2*lst%mynIndex(n) ) = SpectralStateVar(ila,2,k)
2409 SpectralStateRpn(k,2*lst%mymIndex(m) ,2*lst%mynIndex(n)-1) = SpectralStateVar(ila,3,k)
2410 SpectralStateRpn(k,2*lst%mymIndex(m) ,2*lst%mynIndex(n) ) = SpectralStateVar(ila,4,k)
2411 end do
2412 end if
2413 end do
2414 end do
2415 !$OMP END PARALLEL DO
2416
2417 case default
2418 write(*,*)
2419 write(*,*) 'lst_ReshapeTrunc_kij: Unknown Direction', trim(Direction)
2420 call utl_abort('lst_ReshapeTrunc_kij')
2421 end select
2422
2423 end subroutine lst_ReshapeTrunc_kij
2424
2425 !--------------------------------------------------------------------------
2426 ! lst_Laplacian
2427 !--------------------------------------------------------------------------
2428 subroutine lst_Laplacian(lst, GridState, Mode, nk)
2429 implicit none
2430
2431 ! Arguments:
2432 type(struct_lst), intent(in) :: lst
2433 integer, intent(in) :: nk ! Grid point data dimensions
2434 real(8), intent(inout) :: GridState(lst%myLonBeg:lst%myLonEnd,lst%myLatBeg:lst%myLatEnd,nk) ! 3D field in grid point space
2435 character(len=*), intent(in) :: Mode ! Forward or Inverse
2436
2437 ! Locals:
2438 real(8), allocatable :: SpectralStateVar(:,:,:)
2439 real(8), allocatable :: factor(:)
2440 integer :: k, ila, p
2441 character(len=24) :: kind
2442
2443 if (verbose) write(*,*) 'Entering lst_Laplacian'
2444
2445 allocate(SpectralStateVar(lst%nla,lst%nphase,nk))
2446 allocate(factor(lst%nla))
2447
2448 !
2449 !- 1. Set Mode-dependent factors
2450 !
2451 select case (trim(Mode))
2452 case ('Forward')
2453 factor(:) = lst%lapxy(:)
2454 case ('Inverse')
2455 factor(:) = lst%ilapxy(:)
2456 case default
2457 write(*,*)
2458 write(*,*) 'lst_Laplacian: Error: Mode Unknown ', trim(Mode)
2459 call utl_abort('lst_Laplacian')
2460 end select
2461
2462 !
2463 !- 2. Grid Point Space -> Spectral Space
2464 !
2465 kind = 'GridPointToSpectral'
2466 call lst_VarTransform(lst, & ! IN
2467 SpectralStateVar, & ! OUT
2468 GridState, & ! IN
2469 kind, nk) ! IN
2470
2471 !
2472 !- 3. Laplacian (forward or inverse) Transform
2473 !
2474 !$OMP PARALLEL DO PRIVATE (k,ila,p)
2475 do k = 1, nk
2476 do ila = 1, lst%nla
2477 do p = 1, lst%nphase
2478 SpectralStateVar(ila,p,k) = factor(ila) * SpectralStateVar(ila,p,k)
2479 end do
2480 end do
2481 end do
2482 !$OMP END PARALLEL DO
2483
2484 !
2485 !- 4. Spectral Space -> Grid Point Space
2486 !
2487 kind = 'SpectralToGridPoint'
2488 call lst_VarTransform(lst, & ! IN
2489 SpectralStateVar, & ! IN
2490 GridState, & ! OUT
2491 kind, nk) ! IN
2492
2493 deallocate(SpectralStateVar)
2494 deallocate(factor)
2495
2496 end subroutine lst_Laplacian
2497
2498 !--------------------------------------------------------------------------
2499 ! NGFFT
2500 !--------------------------------------------------------------------------
2501 subroutine ngfft(n)
2502 implicit none
2503
2504 ! Arguments:
2505 integer, intent(inout) :: n ! le plus petit entier >= n qui factorise
2506
2507 ! Locals:
2508 integer, parameter :: l = 3
2509 integer :: k(l) , m
2510 data m , k / 8 , 2 , 3 , 5 /
2511 integer :: i, j
2512
2513 if (n <= m) n = m + 1
2514 n = n - 1
25151 n = n + 1
2516 i = n
25172 do j = 1, l
2518 if (mod(i,k(j)) == 0) go to 4
2519 end do
2520 go to 1
25214 i = i/k(j)
2522 if (i /= 1) go to 2
2523
2524 end subroutine ngfft
2525
2526end module lamSpectralTransform_mod
