1module globalSpectralTransform_mod
2 ! MODULE globalSpectralTransform_mod (prefix='gst' category='4. Data Object transformations')
3 !
4 !:Purpose: To perform global spectral transform (spherical harmonic transform
5 ! with grid-point field on a standard global Gaussian grid).
6 !
7 use codePrecision_mod
8 use mpi
9 use midasMpi_mod
10 use MathPhysConstants_mod
11 use earthConstants_mod
12 use utilities_mod
13 implicit none
14 save
15 private
16
17 ! public subroutines
18 public :: gst_setup
19 public :: gst_speree, gst_speree_ad, gst_speree_kij, gst_speree_kij_ad, gst_reespe, gst_reespe_kij
20 public :: gst_spgd, gst_spgda, gst_gdsp, gst_zleginv, gst_zlegdir
21 public :: gst_setID, gst_setDefaultID, gst_setToDefaultID
22 public :: gst_ilaList_mpilocal, gst_ilaList_mpiglobal
23 ! public functions
24 public :: gst_getRmu, gst_getRwt, gst_getnind, gst_getrlati, gst_getr1qm2, gst_getrsqm2
25 public :: gst_getrnnp1, gst_getr1snp1, gst_getzleg, gst_getNla
26
27
28 type :: T_gst
29 real(8),allocatable :: rmu(:)
30 real(8),allocatable :: rwt(:)
31 real(8),allocatable :: rsqm2(:)
32 real(8),allocatable :: r1qm2(:)
33 real(8),allocatable :: rlati(:)
34 integer,allocatable :: nind(:)
35 integer,allocatable :: nindrh(:)
36 real(8),allocatable :: dalp(:,:)
37 real(8),allocatable :: dealp(:,:)
38 real(8),allocatable :: rwocs(:)
39 real(8),allocatable :: r1mu2(:)
40 real(8),allocatable :: rcolat(:)
41 real(8),allocatable :: r1mui(:)
42 real(8),allocatable :: r1mua(:)
43 integer,allocatable :: nclm(:)
44 real(8),allocatable :: r1snp1(:)
45 real(8),allocatable :: rnnp1(:)
46 real(8),allocatable :: zleg(:,:)
47 integer :: ntrunc
48 integer :: nla
49 integer :: nlarh
50 integer :: njlath
51 integer :: ni
52 integer :: nj
53 integer :: nk
54 integer :: myLatBeg, myLatEnd, latPerPE, latPerPEmax
55 integer,allocatable :: allLatBeg(:), allLatEnd(:), allLatPerPE(:)
56 integer :: myLonBeg, myLonEnd, lonPerPE, lonPerPEmax
57 integer,allocatable :: allLonBeg(:), allLonEnd(:), allLonPerPE(:)
58 integer :: myLatHalfBeg, myLatHalfEnd
59 integer :: mymBeg, mymEnd, mymSkip, mymCount, maxmCount
60 integer :: mynBeg, mynEnd, mynSkip, mynCount
61 integer :: myNla, maxMyNla
62 integer,allocatable :: allNla(:)
63 integer,allocatable :: mymIndex(:)
64 integer,pointer :: ilaList(:), allIlaList(:,:)
65 integer,allocatable :: allmBeg(:), allmEnd(:), allmSkip(:)
66 integer,allocatable :: allnBeg(:), allnEnd(:), allnSkip(:)
67 integer :: myLevBeg, myLevEnd, myLevCount, maxMyLevCount
68 integer,allocatable :: allLevBeg(:), allLevEnd(:)
69 integer :: sendType_LevToLon, recvType_LevToLon
70 integer :: sendType_LonToLev, recvType_LonToLev
71 logical :: lonLatDivisible
72 end type T_gst
73
74 integer,parameter :: nMaxGst = 20
75 integer :: gstIdDefault = -1
76 integer :: nGstAlreadyAllocated = 0
77 integer :: gstID = 0
78 type(T_gst) :: gst(nmaxgst)
79
80 integer :: nlatbd = 8
81
82contains
83
84 subroutine GST_setID(gstID_in)
85 implicit none
86
87 ! Arguments:
88 integer, intent(in) :: gstID_in
89
90 gstID = gstID_in
91
92 end subroutine GST_setID
93
94
95 subroutine GST_setDefaultID(gstID_in)
96 implicit none
97
98 ! Arguments:
99 integer, intent(in) :: gstID_in
100
101 gstIDDefault = gstID_in
102
103 end subroutine GST_setDefaultID
104
105
106 subroutine GST_setToDefaultID
107 implicit none
108
109 gstID = gstIdDefault
110
111 end subroutine GST_setToDefaultID
112
113
114 integer function GST_getNla(gstID_opt)
115 implicit none
116
117 ! Arguments:
118 integer, optional, intent(in) :: gstID_opt
119
120 ! Locals:
121 integer :: gstID_l
122
123 if(present(gstID_opt)) then
124 gstID_l = gstID_opt
125 else
126 gstID_l = gstIdDefault
127 endif
128
129 gst_getNla = gst(gstID_l)%myNla
130
131 end function GST_getNla
132
133
134 real(8) function GST_getRmu(latIndex,gstID_opt)
135 implicit none
136
137 ! Arguments:
138 integer, intent(in) :: latIndex
139 integer, optional, intent(in) :: gstID_opt
140
141 ! Locals:
142 integer :: gstID_l
143 integer :: latIndex2
144
145 if(present(gstID_opt)) then
146 gstID_l = gstID_opt
147 else
148 gstID_l = gstIdDefault
149 endif
150
151 ! use a flipped latitude (south to north)
152 latIndex2 = gst(gstID_l)%nj - latIndex + 1
153 gst_getRmu = gst(gstID_l)%rmu(latIndex2)
154
155 end function GST_getRmu
156
157
158 real(8) function GST_getRnnp1(ilaIndex,gstID_opt)
159 implicit none
160
161 ! Arguments:
162 integer, intent(in) :: ilaIndex
163 integer, optional, intent(in) :: gstID_opt
164
165 ! Locals:
166 integer :: gstID_l
167
168 if(present(gstID_opt)) then
169 gstID_l = gstID_opt
170 else
171 gstID_l = gstIdDefault
172 endif
173
174 gst_getRnnp1 = gst(gstID_l)%rnnp1(ilaIndex)
175
176 end function GST_getRnnp1
177
178
179 real(8) function GST_getR1snp1(ilaIndex,gstID_opt)
180 implicit none
181
182 ! Arguments:
183 integer, intent(in) :: ilaIndex
184 integer, optional, intent(in) :: gstID_opt
185
186 ! Locals:
187 integer :: gstID_l
188
189 if(present(gstID_opt)) then
190 gstID_l = gstID_opt
191 else
192 gstID_l = gstIdDefault
193 endif
194
195 gst_getR1snp1 = gst(gstID_l)%r1snp1(ilaIndex)
196
197 end function GST_getR1snp1
198
199
200 real(8) function GST_getRwt(latIndex,gstID_opt)
201 implicit none
202
203 ! Arguments:
204 integer, intent(in) :: latIndex
205 integer, optional, intent(in) :: gstID_opt
206
207 ! Locals:
208 integer :: gstID_l
209 integer :: latIndex2
210
211 if(present(gstID_opt)) then
212 gstID_l = gstID_opt
213 else
214 gstID_l = gstIdDefault
215 endif
216
217 ! use a flipped latitude (south to north)
218 latIndex2 = gst(gstID_l)%nj - latIndex + 1
219 gst_getRwt = gst(gstID_l)%rwt(latIndex2)
220
221 end function GST_getRwt
222
223
224 integer function GST_getNind(mIndex,gstID_opt)
225 implicit none
226
227 ! Arguments:
228 integer, intent(in) :: mIndex
229 integer, optional, intent(in) :: gstID_opt
230
231 ! Locals:
232 integer :: gstID_l
233
234 if(present(gstID_opt)) then
235 gstID_l = gstID_opt
236 else
237 gstID_l = gstIdDefault
238 endif
239
240 gst_getNind = gst(gstID_l)%nind(mIndex)
241
242 end function GST_getNind
243
244
245 real(8) function GST_getRlati(latIndex,gstID_opt)
246 implicit none
247
248 ! Arguments:
249 integer, intent(in) :: latIndex
250 integer, optional, intent(in) :: gstID_opt
251
252 ! Locals:
253 integer :: gstID_l
254 integer :: latIndex2
255
256 if(present(gstID_opt)) then
257 gstID_l = gstID_opt
258 else
259 gstID_l = gstIdDefault
260 endif
261
262 ! use a flipped latitude (south to north)
263 latIndex2 = gst(gstID_l)%nj - latIndex + 1
264 gst_getRlati = gst(gstID_l)%rlati(latIndex2)
265
266 end function GST_getRlati
267
268
269 real(8) function GST_getR1qm2(latIndex,gstID_opt)
270 implicit none
271
272 ! Arguments:
273 integer, intent(in) :: latIndex
274 integer, optional, intent(in) :: gstID_opt
275
276 ! Locals:
277 integer :: gstID_l
278 integer :: latIndex2
279
280 if(present(gstID_opt)) then
281 gstID_l = gstID_opt
282 else
283 gstID_l = gstIdDefault
284 endif
285
286 ! use a flipped latitude (south to north)
287 latIndex2 = gst(gstID_l)%nj - latIndex + 1
288 gst_getR1qm2 = gst(gstID_l)%r1qm2(latIndex2)
289
290 end function GST_getR1qm2
291
292
293 real(8) function GST_getRsqm2(latIndex,gstID_opt)
294 implicit none
295
296 ! Arguments:
297 integer, intent(in) :: latIndex
298 integer, optional, intent(in) :: gstID_opt
299
300 ! Locals:
301 integer :: gstID_l
302 integer :: latIndex2
303
304 if(present(gstID_opt)) then
305 gstID_l = gstID_opt
306 else
307 gstID_l = gstIdDefault
308 endif
309
310 ! use a flipped latitude (south to north)
311 latIndex2 = gst(gstID_l)%nj - latIndex + 1
312 gst_getRsqm2 = gst(gstID_l)%rsqm2(latIndex2)
313
314 end function GST_getRsqm2
315
316
317 real(8) function GST_getzleg(legendreIndex,latIndex,gstID_in)
318 !
319 !:Purpose: To pass on Legendre polynomial element
320 !
321 implicit none
322
323 ! Arguments:
324 integer :: legendreIndex
325 integer :: latIndex
326 integer :: gstID_in
327
328 gst_getzleg=gst(gstID_in)%zleg(legendreIndex,latIndex)
329
330 end function GST_getzleg
331
332
333 subroutine GST_ilaList_mpiglobal(ilaList,myNla,maxMyNla,gstID_in,mymBeg,&
334 mymEnd,mymSkip,mynBeg,mynEnd,mynSkip)
335 !
336 !:Purpose: To produce an array to convert an mpilocal "ila" into an
337 ! mpiglobal "ila"
338 implicit none
339
340 ! Arguments:
341 integer, pointer, intent(out) :: ilaList(:)
342 integer, intent(out) :: myNla
343 integer, intent(in) :: maxMyNla
344 integer, intent(in) :: gstID_in
345 integer, intent(in) :: mymBeg
346 integer, intent(in) :: mymEnd
347 integer, intent(in) :: mymSkip
348 integer, intent(in) :: mynBeg
349 integer, intent(in) :: mynEnd
350 integer, intent(in) :: mynSkip
351
352 ! Locals:
353 integer :: jm, jn, ierr
354
355 ! compute mpilocal value of nla
356 myNla = 0
357 do jm = mymBeg, mymEnd, mymSkip
358 do jn = mynBeg, mynEnd, mynSkip
359 if(jm.le.jn) then
360 myNla = myNla + 1
361 endif
362 enddo
363 enddo
364
365 ! determine maximum value of myNla over all processors (used for dimensioning)
366 call rpn_comm_allreduce(myNla,maxMyNla,1,'MPI_INTEGER','MPI_MAX','GRID',ierr)
367
368 allocate(ilaList(maxMyNla))
369 ilaList(:) = 0
370
371 myNla = 0
372 do jm = mymBeg, mymEnd, mymSkip
373 do jn = mynBeg, mynEnd, mynSkip
374 if(jm.le.jn) then
375 myNla = myNla + 1
376 ilaList(myNla) = gst_getNind(jm,gstID_in) + jn - jm
377 endif
378 enddo
379 enddo
380
381 end subroutine GST_ilaList_mpiglobal
382
383
384 subroutine GST_ilaList_mpilocal(ilaList,gstID_in,mymBeg,mymEnd,mymSkip,&
385 mynBeg,mynEnd,mynSkip)
386 !
387 !:Purpose: To produce an array to convert an mpiglobal "ila" into an
388 ! mpilocal "ila"
389 !
390 implicit none
391
392 ! Arguments:
393 integer, pointer, intent(out) :: ilaList(:)
394 integer, intent(in) :: gstID_in
395 integer, intent(in) :: mymBeg
396 integer, intent(in) :: mymEnd
397 integer, intent(in) :: mymSkip
398 integer, intent(in) :: mynBeg
399 integer, intent(in) :: mynEnd
400 integer, intent(in) :: mynSkip
401
402 ! Locals:
403 integer :: jm, jn, myNla
404
405 ! assume mpiglobal value of nla already set in gst structure
406 allocate(ilaList(gst(gstID_in)%nla))
407 ilaList(:) = 0
408
409 myNla = 0
410 do jm = mymBeg, mymEnd, mymSkip
411 do jn = mynBeg, mynEnd, mynSkip
412 if(jm.le.jn) then
413 myNla = myNla + 1
414 ilaList(gst_getNind(jm,gstID_in) + jn - jm) = myNla
415 endif
416 enddo
417 enddo
418
419 end subroutine GST_ilaList_mpilocal
420
421
422 integer function gst_setup(ni_in,nj_in,ntrunc_in,maxlevels_in)
423 implicit none
424
425 ! Arguments:
426 integer, intent(in) :: ni_in
427 integer, intent(in) :: nj_in
428 integer, intent(in) :: ntrunc_in
429 integer, intent(in) :: maxlevels_in
430
431 ! Locals:
432 integer :: jn, jm, ila, ierr
433 integer :: latPerPE, latPerPEmax, myLatBeg, myLatEnd, myLatHalfBeg, myLatHalfEnd
434 integer :: lonPerPE, lonPerPEmax, myLonBeg, myLonEnd
435 integer :: myLevBeg, myLevEnd, myLevCount
436 integer :: mymBeg, mymEnd, mymSkip, mymCount
437 integer :: mynBeg, mynEnd, mynSkip, mynCount
438 real(8) :: znnp1, z1snp1
439 integer(kind=MPI_ADDRESS_KIND) :: lowerBound, extent
440 integer :: realSize, sendType, recvType
441 logical :: divisibleLon, divisibleLat
442
443 if(nGstAlreadyAllocated.eq.nMaxGst) then
444 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: The maxmimum number of spectral transform have already been allocated! ', nMaxGst
445 call utl_abort('gst_setup')
446 endif
447
448 nGstAlreadyAllocated = nGstAlreadyAllocated+1
449 gstID = nGstAlreadyAllocated
450 call gst_setDefaultID(gstID)
451 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: Now setting up spectral transform #', gstID
452 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: following *kind* used for internal reals : ', pre_specTransReal
453
454 gst(gstID)%ni = ni_in
455 gst(gstID)%nj = nj_in
456 gst(gstID)%njlath = (gst(gstID)%nj + 1)/2
457
458 if (ntrunc_in == -1) then ! no truncation case
459 gst(gstID)%ntrunc = gst(gstID)%ni / 2
460 else
461 gst(gstID)%ntrunc = ntrunc_in
462 end if
463 gst(gstID)%nla = (gst(gstID)%ntrunc + 1) * (gst(gstID)%ntrunc +2)/2
464 gst(gstID)%nlarh = (gst(gstID)%ntrunc+1) * (gst(gstID)%ntrunc+1)
465 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: ntrunc=', gst(gstID)%ntrunc
466
467 call mmpi_setup_latbands(gst(gstID)%nj, &
468 latPerPE, latPerPEmax, myLatBeg, myLatEnd, myLatHalfBeg, myLatHalfEnd, divisible_opt=divisibleLat)
469 call mmpi_setup_lonbands(gst(gstID)%ni, &
470 lonPerPE, lonPerPEmax, myLonBeg, myLonEnd, divisible_opt= divisibleLon)
471
472 gst(gstID)%lonLatDivisible = (divisibleLon .and. divisibleLat)
473 if( mmpi_myid == 0 ) write(*,*) 'gst_setup: lonLatDivisible = ', gst(gstID)%lonLatDivisible
474
475 gst(gstID)%nk = maxlevels_in
476 ! 2D MPI decomposition: split levels across npex
477 call mmpi_setup_levels(maxlevels_in, myLevBeg, myLevEnd, myLevCount)
478 write(*,*) 'gst_setup: myLevBeg,End,Count=', myLevBeg, myLevEnd, myLevCount
479
480 !!! Distribution of lon/lat tiles (gridpoint space) and n/m (spectral space)
481 ! range of lons handled by this processor
482 gst(gstID)%myLonBeg = myLonBeg
483 gst(gstID)%myLonEnd = myLonEnd
484 gst(gstID)%lonPerPE = lonPerPE
485 gst(gstID)%lonPerPEmax = lonPerPEmax
486 ! range of lats handled by this processor
487 gst(gstID)%myLatBeg = myLatBeg
488 gst(gstID)%myLatEnd = myLatEnd
489 gst(gstID)%myLatHalfBeg = myLatHalfBeg
490 gst(gstID)%myLatHalfEnd = myLatHalfEnd
491 gst(gstID)%latPerPE = latPerPE
492 gst(gstID)%latPerPEmax = latPerPEmax
493 ! range of n handled by this processor
494 call mmpi_setup_n(gst(gstID)%ntrunc,mynBeg,mynEnd,mynSkip,mynCount)
495 gst(gstID)%mynBeg = mynBeg
496 gst(gstID)%mynEnd = mynEnd
497 gst(gstID)%mynSkip = mynSkip
498 gst(gstID)%mynCount = mynCount
499 ! range of m handled by this processor
500 call mmpi_setup_m(gst(gstID)%ntrunc,mymBeg,mymEnd,mymSkip,mymCount)
501 gst(gstID)%mymBeg = mymBeg
502 gst(gstID)%mymEnd = mymEnd
503 gst(gstID)%mymSkip = mymSkip
504 gst(gstID)%mymCount = mymCount
505 call rpn_comm_allreduce(gst(gstID)%mymCount,gst(gstID)%maxmCount, &
506 1,'MPI_INTEGER','MPI_MAX','GRID',ierr)
507 ! range of levels handled by this processor when in spectral space
508 gst(gstID)%myLevBeg = myLevBeg
509 gst(gstID)%myLevEnd = myLevEnd
510 gst(gstID)%myLevCount = myLevCount
511 call rpn_comm_allreduce(gst(gstID)%myLevCount,gst(gstID)%maxMyLevCount, &
512 1,'MPI_INTEGER','MPI_MAX','GRID',ierr)
513
514 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allocating comleg...'
515 call allocate_comleg
516 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: calling suleg...'
517 call suleg(.false.)
518 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: calling sualp...'
519 call sualp
520
521 ! compute the list of spectral coefficient indices when distributed by n and m
522 call gst_ilaList_mpiglobal(gst(gstID)%ilaList,gst(gstID)%myNla,gst(gstID)%maxMyNla,gstID, &
523 gst(gstID)%mymBeg,gst(gstID)%mymEnd,gst(gstID)%mymSkip, &
524 gst(gstID)%mynBeg,gst(gstID)%mynEnd,gst(gstID)%mynSkip)
525 allocate(gst(gstID)%allNla(mmpi_npex))
526 call rpn_comm_allgather(gst(gstID)%myNla,1,'mpi_integer', &
527 gst(gstID)%allNla,1,'mpi_integer','EW',ierr)
528 allocate(gst(gstID)%allIlaList(gst(gstID)%maxMyNla,mmpi_npex))
529 call rpn_comm_allgather(gst(gstID)%ilaList,gst(gstID)%maxMyNla,'mpi_integer', &
530 gst(gstID)%allIlaList,gst(gstID)%maxMyNla,'mpi_integer','EW',ierr)
531
532 allocate(gst(gstID)%allLonBeg(mmpi_npex))
533 CALL rpn_comm_allgather(gst(gstID)%myLonBeg,1,'mpi_integer', &
534 gst(gstID)%allLonBeg,1,'mpi_integer','EW',ierr)
535 allocate(gst(gstID)%allLonEnd(mmpi_npex))
536 CALL rpn_comm_allgather(gst(gstID)%myLonEnd,1,'mpi_integer', &
537 gst(gstID)%allLonEnd,1,'mpi_integer','EW',ierr)
538 allocate(gst(gstID)%allLonPerPE(mmpi_npex))
539 CALL rpn_comm_allgather(gst(gstID)%lonPerPE,1,'mpi_integer', &
540 gst(gstID)%allLonPerPE,1,'mpi_integer','EW',ierr)
541
542 allocate(gst(gstID)%allLatBeg(mmpi_npey))
543 CALL rpn_comm_allgather(gst(gstID)%myLatBeg,1,'mpi_integer', &
544 gst(gstID)%allLatBeg,1,'mpi_integer','NS',ierr)
545 allocate(gst(gstID)%allLatEnd(mmpi_npey))
546 CALL rpn_comm_allgather(gst(gstID)%myLatEnd,1,'mpi_integer', &
547 gst(gstID)%allLatEnd,1,'mpi_integer','NS',ierr)
548 allocate(gst(gstID)%allLatPerPE(mmpi_npey))
549 CALL rpn_comm_allgather(gst(gstID)%latPerPE,1,'mpi_integer', &
550 gst(gstID)%allLatPerPE,1,'mpi_integer','NS',ierr)
551
552 allocate(gst(gstID)%mymIndex(gst(gstID)%mymBeg:gst(gstID)%mymEnd))
553 gst(gstID)%mymIndex(:) = 0
554 do jm = gst(gstID)%mymBeg,gst(gstID)%mymEnd,gst(gstID)%mymSkip
555 if(jm.eq.gst(gstID)%mymBeg) then
556 gst(gstID)%mymIndex(jm) = 1
557 else
558 gst(gstID)%mymIndex(jm) = gst(gstID)%mymIndex(jm-gst(gstID)%mymSkip) + 1
559 endif
560 write(*,*) 'gst_setup: mymIndex(',jm,')=',gst(gstID)%mymIndex(jm)
561 enddo
562
563 allocate(gst(gstID)%allnBeg(mmpi_npex))
564 CALL rpn_comm_allgather(gst(gstID)%mynBeg,1,'mpi_integer', &
565 gst(gstID)%allnBeg,1,'mpi_integer','EW',ierr)
566 allocate(gst(gstID)%allnEnd(mmpi_npex))
567 CALL rpn_comm_allgather(gst(gstID)%mynEnd,1,'mpi_integer', &
568 gst(gstID)%allnEnd,1,'mpi_integer','EW',ierr)
569 allocate(gst(gstID)%allnSkip(mmpi_npex))
570 CALL rpn_comm_allgather(gst(gstID)%mynSkip,1,'mpi_integer', &
571 gst(gstID)%allnSkip,1,'mpi_integer','EW',ierr)
572
573 allocate(gst(gstID)%allmBeg(mmpi_npey))
574 CALL rpn_comm_allgather(gst(gstID)%mymBeg,1,'mpi_integer', &
575 gst(gstID)%allmBeg,1,'mpi_integer','NS',ierr)
576 allocate(gst(gstID)%allmEnd(mmpi_npey))
577 CALL rpn_comm_allgather(gst(gstID)%mymEnd,1,'mpi_integer', &
578 gst(gstID)%allmEnd,1,'mpi_integer','NS',ierr)
579 allocate(gst(gstID)%allmSkip(mmpi_npey))
580 CALL rpn_comm_allgather(gst(gstID)%mymSkip,1,'mpi_integer', &
581 gst(gstID)%allmSkip,1,'mpi_integer','NS',ierr)
582
583 allocate(gst(gstID)%allLevBeg(mmpi_npex))
584 CALL rpn_comm_allgather(gst(gstID)%myLevBeg,1,'mpi_integer', &
585 gst(gstID)%allLevBeg,1,'mpi_integer','EW',ierr)
586 allocate(gst(gstID)%allLevEnd(mmpi_npex))
587 CALL rpn_comm_allgather(gst(gstID)%myLevEnd,1,'mpi_integer', &
588 gst(gstID)%allLevEnd,1,'mpi_integer','EW',ierr)
589
590 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allLonBeg=',gst(gstID)%allLonBeg
591 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allLonEnd=',gst(gstID)%allLonEnd
592 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allLatBeg=',gst(gstID)%allLatBeg
593 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allLatEnd=',gst(gstID)%allLatEnd
594 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allnBeg=',gst(gstID)%allnBeg
595 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allnEnd=',gst(gstID)%allnEnd
596 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allnSkip=',gst(gstID)%allnSkip
597 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allmBeg=',gst(gstID)%allmBeg
598 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allmEnd=',gst(gstID)%allmEnd
599 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allmSkip=',gst(gstID)%allmSkip
600 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allLevBeg=',gst(gstID)%allLevBeg
601 if(mmpi_myid.eq.0) write(*,*) 'gst_setup: allLevEnd=',gst(gstID)%allLevEnd
602 write(*,*) 'gst_setup: mymCount=',gst(gstID)%mymCount
603 write(*,*) 'gst_setup: maxmCount=',gst(gstID)%maxmCount
604 write(*,*) 'gst_setup: myNla=',gst(gstID)%myNla
605 write(*,*) 'gst_setup: maxMyNla=',gst(gstID)%maxMyNla
606
607 allocate(gst(gstID)%r1snp1(gst(gstID)%nla))
608 allocate(gst(gstID)%rnnp1(gst(gstID)%nla))
609 gst(gstID)%r1snp1(1) = 0.d0
610 gst(gstID)%rnnp1(1) = 0.d0
611 do jn = 1, gst(gstID)%ntrunc
612 znnp1 = -real(jn,8)*real(jn+1,8)
613 z1snp1 = 1.d0/znnp1
614 do jm = 0, jn
615 ila = gst(gstID)%nind(jm) + jn - jm
616 gst(gstID)%r1snp1(ila) = z1snp1
617 gst(gstID)%rnnp1(ila) = znnp1
618 enddo
619 enddo
620
621 call gst_zlegpol(gstID)
622
623 ! setup mpi derived types used in transposes (only used when grid is divisible)
624 ! ... mpi_type_vector(count, blocklength, stride, ...)
625 ! ... mpi_type_create_resized(oldtype, lowerbound, extent(in bytes), newtype, ierr)
626
627 call mpi_type_size(pre_specTransMpiType, realSize, ierr)
628 lowerBound = 0
629
630 ! create the send type for LevToLon
631 extent = gst(gstID)%maxMyLevCount * gst(gstID)%lonPerPE * realSize
632 call mpi_type_vector(gst(gstID)%latPerPE, gst(gstID)%maxMyLevCount * gst(gstID)%lonPerPE, &
633 gst(gstID)%maxMyLevCount * gst(gstID)%ni, pre_specTransMpiType, sendtype, ierr)
634 call mpi_type_create_resized(sendtype, lowerBound , extent, gst(gstID)%sendType_LevToLon, ierr);
635 call mpi_type_commit(gst(gstID)%sendType_LevToLon,ierr)
636
637 ! create the receive type for LevToLon
638 extent = gst(gstID)%maxMyLevCount * realSize
639 call mpi_type_vector(gst(gstID)%lonPerPE * gst(gstID)%latPerPE , gst(gstID)%maxMyLevCount, &
640 gst(gstID)%nk, pre_specTransMpiType, recvtype, ierr);
641 call mpi_type_create_resized(recvtype, lowerBound, extent, gst(gstID)%recvType_LevToLon, ierr);
642 call mpi_type_commit(gst(gstID)%recvType_LevToLon, ierr);
643
644 ! create the send type for LonToLev
645 extent = gst(gstID)%maxMyLevCount * realSize
646 call mpi_type_vector(gst(gstID)%lonPerPE * gst(gstID)%latPerPE , gst(gstID)%maxMyLevCount, &
647 gst(gstID)%nk, pre_specTransMpiType, sendtype, ierr);
648 call mpi_type_create_resized(sendtype, lowerBound, extent, gst(gstID)%sendType_LonToLev, ierr);
649 call mpi_type_commit(gst(gstID)%sendType_LonToLev, ierr);
650
651 ! create the recv type for LonToLev
652 extent = gst(gstID)%maxMyLevCount * gst(gstID)%lonPerPE * realSize
653 call mpi_type_vector(gst(gstID)%latPerPE, gst(gstID)%maxMyLevCount * gst(gstID)%lonPerPE, &
654 gst(gstID)%maxMyLevCount * gst(gstID)%ni, pre_specTransMpiType, recvtype, ierr)
655 call mpi_type_create_resized(recvtype, lowerBound , extent, gst(gstID)%recvType_LonToLev, ierr);
656 call mpi_type_commit(gst(gstID)%recvType_LonToLev,ierr)
657
658 gst_setup = gstID
659
660 end function GST_SETUP
661
662!-------------------------------------------------------------------------------
663! Data transposes with respect to 1 of 2 dimensions (i.e. NS or EW)
664!-------------------------------------------------------------------------------
665
666 subroutine transpose2d_NtoLev(psp_in,psp_out)
667 implicit none
668
669 ! Arguments:
670 real(8), intent(in) :: psp_in (gst(gstID)%myNla, 2, gst(gstID)%nk)
671 real(8), intent(out) :: psp_out(gst(gstID)%nla, 2, &
672 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
673
674 ! Locals:
675 real(pre_specTransReal) :: sp_send(gst(gstID)%maxMyNla, 2, &
676 gst(gstID)%maxMyLevCount, mmpi_npex)
677 real(pre_specTransReal) :: sp_recv(gst(gstID)%maxMyNla, 2, &
678 gst(gstID)%maxMyLevCount, mmpi_npex)
679 integer :: yourid,ila,icount,nsize,ierr,jlev,jlev2
680
681 call utl_tmg_start(152,'low-level--gst_barr')
682 if(mmpi_doBarrier) call rpn_comm_barrier('EW',ierr)
683 call utl_tmg_stop(152)
684
685 call utl_tmg_start(153,'low-level--gst_transpose_NtoLEV')
686
687 !$OMP PARALLEL DO PRIVATE(yourid,jlev,jlev2,icount)
688 do yourid = 0, (mmpi_npex-1)
689 do jlev = gst(gstID)%allLevBeg(yourid+1), gst(gstID)%allLevEnd(yourid+1)
690 jlev2 = jlev - gst(gstID)%allLevBeg(yourid+1) + 1
691 do icount = 1, gst(gstID)%myNla
692 sp_send(icount,1,jlev2,yourid+1) = psp_in(icount,1,jlev)
693 sp_send(icount,2,jlev2,yourid+1) = psp_in(icount,2,jlev)
694 enddo
695 enddo
696 enddo
697 !$OMP END PARALLEL DO
698
699 nsize = gst(gstID)%maxMyNla * 2 * gst(gstID)%maxMyLevCount
700 if(mmpi_npex.gt.1) then
701 call rpn_comm_alltoall(sp_send, nsize, pre_specTransMpiReal, &
702 sp_recv, nsize, pre_specTransMpiReal, 'EW', ierr)
703 else
704 sp_recv(:,:,:,1) = sp_send(:,:,:,1)
705 endif
706
707 !$OMP PARALLEL DO PRIVATE(yourid,jlev,jlev2,icount,ila)
708 do yourid = 0, (mmpi_npex-1)
709 do jlev = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
710 jlev2 = jlev - gst(gstID)%myLevBeg + 1
711 do icount = 1, gst(gstID)%allNla(yourid+1)
712 ila = gst(gstID)%allIlaList(icount,yourid+1)
713 psp_out(ila,1,jlev) = sp_recv(icount,1,jlev2,yourid+1)
714 psp_out(ila,2,jlev) = sp_recv(icount,2,jlev2,yourid+1)
715 enddo
716 enddo
717 enddo
718 !$OMP END PARALLEL DO
719
720 call utl_tmg_stop(153)
721
722 end subroutine transpose2d_NtoLev
723
724
725 subroutine transpose2d_LevtoN(psp_in,psp_out)
726 implicit none
727
728 ! Arguments:
729 real(8), intent(in) :: psp_in (gst(gstID)%nla, 2, &
730 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
731 real(8), intent(out) :: psp_out(gst(gstID)%myNla, 2, gst(gstID)%nk)
732
733 ! Locals:
734 real(pre_specTransReal) :: sp_send(gst(gstID)%maxMyNla, 2, &
735 gst(gstID)%maxMyLevCount, mmpi_npex)
736 real(pre_specTransReal) :: sp_recv(gst(gstID)%maxMyNla, 2, &
737 gst(gstID)%maxMyLevCount, mmpi_npex)
738 integer :: yourid,ila,icount,nsize,ierr,jlev,jlev2
739
740 call utl_tmg_start(152,'low-level--gst_barr')
741 if(mmpi_doBarrier) call rpn_comm_barrier('EW',ierr)
742 call utl_tmg_stop(152)
743
744 call utl_tmg_start(153,'low-level--gst_transpose_NtoLEV')
745
746 !$OMP PARALLEL DO PRIVATE(yourid,jlev,jlev2,icount,ila)
747 do yourid = 0, (mmpi_npex-1)
748 do jlev = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
749 jlev2 = jlev - gst(gstID)%myLevBeg + 1
750 sp_send(:,:,jlev2,yourid+1) = 0.0d0
751 do icount = 1, gst(gstID)%allNla(yourid+1)
752 ila = gst(gstID)%allIlaList(icount,yourid+1)
753 sp_send(icount,1,jlev2,yourid+1) = psp_in(ila,1,jlev)
754 sp_send(icount,2,jlev2,yourid+1) = psp_in(ila,2,jlev)
755 enddo
756 enddo
757 enddo
758 !$OMP END PARALLEL DO
759
760 nsize = gst(gstID)%maxMyNla * 2 * gst(gstID)%maxMyLevCount
761 if(mmpi_npex.gt.1) then
762 call rpn_comm_alltoall(sp_send, nsize, pre_specTransMpiReal, &
763 sp_recv, nsize, pre_specTransMpiReal, 'EW', ierr)
764 else
765 sp_recv(:,:,:,1) = sp_send(:,:,:,1)
766 endif
767
768 !$OMP PARALLEL DO PRIVATE(yourid,jlev,jlev2,icount)
769 do yourid = 0, (mmpi_npex-1)
770 do jlev = gst(gstID)%allLevBeg(yourid+1), gst(gstID)%allLevEnd(yourid+1)
771 jlev2 = jlev - gst(gstID)%allLevBeg(yourid+1) + 1
772 do icount = 1, gst(gstID)%myNla
773 psp_out(icount,1,jlev) = sp_recv(icount,1,jlev2,yourid+1)
774 psp_out(icount,2,jlev) = sp_recv(icount,2,jlev2,yourid+1)
775 enddo
776 enddo
777 enddo
778 !$OMP END PARALLEL DO
779
780 call utl_tmg_stop(153)
781
782 end subroutine transpose2d_LevtoN
783
784
785 subroutine transpose2d_MtoLat(pgd_in,pgd_out)
786 implicit none
787
788 ! Arguments:
789 real(8), intent(in) :: pgd_in(2*gst(gstID)%maxmCount, gst(gstID)%nj, &
790 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
791 real(8), intent(out) :: pgd_out(gst(gstID)%ni, gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, &
792 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
793
794 ! Locals:
795 real(pre_specTransReal) :: gd_send(gst(gstID)%maxmCount, 2, gst(gstID)%latPerPEmax, &
796 gst(gstID)%maxMyLevCount, mmpi_npey)
797 real(pre_specTransReal) :: gd_recv(gst(gstID)%maxmCount, 2, gst(gstID)%latPerPEmax, &
798 gst(gstID)%maxMyLevCount, mmpi_npey)
799 integer :: yourid,jm,jm2,icount,nsize,ierr,jlev,jlev2,jlat,jlat2
800
801 call utl_tmg_start(152,'low-level--gst_barr')
802 if(mmpi_doBarrier) call rpn_comm_barrier('NS',ierr)
803 call utl_tmg_stop(152)
804
805 call utl_tmg_start(154,'low-level--gst_transpose_MtoLAT')
806
807 !$OMP PARALLEL DO PRIVATE(yourid,jlat,jlat2,jlev,jlev2,icount,jm,jm2)
808 do yourid = 0, (mmpi_npey-1)
809 do jlev = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
810 jlev2 = jlev - gst(gstID)%myLevBeg + 1
811 gd_send(:, :, :, jlev2, yourid+1) = 0.0d0
812 do jlat = gst(gstID)%allLatBeg(yourid+1), gst(gstID)%allLatEnd(yourid+1)
813 jlat2 = jlat - gst(gstID)%allLatBeg(yourid+1) + 1
814 icount = 0
815 do jm = gst(gstID)%mymBeg, gst(gstID)%mymEnd, gst(gstID)%mymSkip
816 jm2 = 2*gst(gstID)%mymIndex(jm)
817 icount = icount + 1
818 gd_send(icount, 1, jlat2, jlev2, yourid+1) = pgd_in(jm2-1, jlat, jlev)
819 gd_send(icount, 2, jlat2, jlev2, yourid+1) = pgd_in(jm2 , jlat, jlev)
820 enddo
821 enddo
822 enddo
823 enddo
824 !$OMP END PARALLEL DO
825
826 nsize = gst(gstID)%maxmCount * 2 * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPEmax
827 if(mmpi_npey.gt.1) then
828 call rpn_comm_alltoall(gd_send, nsize, pre_specTransMpiReal, &
829 gd_recv, nsize, pre_specTransMpiReal, 'NS', ierr)
830 else
831 gd_recv(:,:,:,:,1) = gd_send(:,:,:,:,1)
832 endif
833
834 !$OMP PARALLEL DO PRIVATE(yourid,jlat,jlat2,jlev,jlev2,icount,jm,jm2)
835 do yourid = 0, (mmpi_npey-1)
836 do jlev = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
837 jlev2 = jlev - gst(gstID)%myLevBeg + 1
838 do jlat = gst(gstID)%myLatBeg, gst(gstID)%myLatEnd
839 jlat2 = jlat - gst(gstID)%myLatBeg + 1
840 icount = 0
841 do jm = gst(gstID)%allmBeg(yourid+1), gst(gstID)%allmEnd(yourid+1), gst(gstID)%allmSkip(yourid+1)
842 jm2 = 2*jm+1
843 icount = icount + 1
844 pgd_out(jm2 , jlat, jlev) = gd_recv(icount, 1, jlat2, jlev2, yourid+1)
845 pgd_out(jm2+1, jlat, jlev) = gd_recv(icount, 2, jlat2, jlev2, yourid+1)
846 enddo
847 enddo
848 enddo
849 enddo
850 !$OMP END PARALLEL DO
851
852 call utl_tmg_stop(154)
853
854 end subroutine transpose2d_MtoLat
855
856
857 subroutine transpose2d_MtoLat_kij(pgd_in,pgd_out)
858 implicit none
859
860 ! Arguments:
861 real(8), intent(in) :: pgd_in(gst(gstID)%maxMyLevCount, 2*gst(gstID)%maxmCount, &
862 gst(gstID)%nj)
863 real(8), intent(out) :: pgd_out(gst(gstID)%maxMyLevCount, gst(gstID)%ni, &
864 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
865
866 ! Locals:
867 real(pre_specTransReal), allocatable, save :: gd_send(:,:,:,:,:)
868 real(pre_specTransReal), allocatable, save :: gd_recv(:,:,:,:,:)
869 integer :: yourid,jm,jm2,icount,nsize,ierr,jlev,jlat,jlat2
870
871 call utl_reAllocate(gd_send, gst(gstID)%maxMyLevCount, gst(gstID)%maxmCount, 2, gst(gstID)%latPerPEmax, mmpi_npey)
872 call utl_reAllocate(gd_recv, gst(gstID)%maxMyLevCount, gst(gstID)%maxmCount, 2, gst(gstID)%latPerPEmax, mmpi_npey)
873
874 call utl_tmg_start(152,'low-level--gst_barr')
875 if(mmpi_doBarrier) call rpn_comm_barrier('NS',ierr)
876 call utl_tmg_stop(152)
877
878 call utl_tmg_start(154,'low-level--gst_transpose_MtoLAT')
879
880 !$OMP PARALLEL DO PRIVATE(yourid,jlat,jlat2,jlev,icount,jm,jm2)
881 do yourid = 0, (mmpi_npey-1)
882 gd_send(:, :, :, :, yourid+1) = 0.0d0
883 do jlat = gst(gstID)%allLatBeg(yourid+1), gst(gstID)%allLatEnd(yourid+1)
884 jlat2 = jlat - gst(gstID)%allLatBeg(yourid+1) + 1
885 icount = 0
886 do jm = gst(gstID)%mymBeg, gst(gstID)%mymEnd, gst(gstID)%mymSkip
887 jm2 = 2*gst(gstID)%mymIndex(jm)
888 icount = icount + 1
889 do jlev = 1, gst(gstID)%maxMyLevCount
890 gd_send(jlev, icount, 1, jlat2, yourid+1) = pgd_in(jlev, jm2-1, jlat)
891 gd_send(jlev, icount, 2, jlat2, yourid+1) = pgd_in(jlev, jm2 , jlat)
892 enddo
893 enddo
894 enddo
895 enddo
896 !$OMP END PARALLEL DO
897
898 nsize = gst(gstID)%maxmCount * 2 * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPEmax
899 if(mmpi_npey.gt.1) then
900 call rpn_comm_alltoall(gd_send, nsize, pre_specTransMpiReal, &
901 gd_recv, nsize, pre_specTransMpiReal, 'NS', ierr)
902 else
903 gd_recv(:,:,:,:,1) = gd_send(:,:,:,:,1)
904 endif
905
906 !$OMP PARALLEL DO PRIVATE(yourid,jlat,jlat2,jlev,icount,jm,jm2)
907 do yourid = 0, (mmpi_npey-1)
908 do jlat = gst(gstID)%myLatBeg, gst(gstID)%myLatEnd
909 jlat2 = jlat - gst(gstID)%myLatBeg + 1
910 icount = 0
911 do jm = gst(gstID)%allmBeg(yourid+1), gst(gstID)%allmEnd(yourid+1), gst(gstID)%allmSkip(yourid+1)
912 jm2 = 2*jm+1
913 icount = icount + 1
914 do jlev = 1, gst(gstID)%maxMyLevCount
915 pgd_out(jlev, jm2 , jlat) = gd_recv(jlev, icount, 1, jlat2, yourid+1)
916 pgd_out(jlev, jm2+1, jlat) = gd_recv(jlev, icount, 2, jlat2, yourid+1)
917 enddo
918 enddo
919 enddo
920 enddo
921 !$OMP END PARALLEL DO
922
923 call utl_tmg_stop(154)
924
925 end subroutine transpose2d_MtoLat_kij
926
927
928 subroutine transpose2d_LattoM(pgd_in,pgd_out)
929 implicit none
930
931 ! Arguments:
932 real(8), intent(in) :: pgd_in(gst(gstID)%ni, gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, &
933 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
934 real(8), intent(out) :: pgd_out(2*gst(gstID)%maxmCount, gst(gstID)%nj, &
935 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
936
937 ! Locals:
938 real(pre_specTransReal) :: gd_send(gst(gstID)%maxmCount, 2, gst(gstID)%latPerPEmax, &
939 gst(gstID)%maxMyLevCount, mmpi_npey)
940 real(pre_specTransReal) :: gd_recv(gst(gstID)%maxmCount, 2, gst(gstID)%latPerPEmax, &
941 gst(gstID)%maxMyLevCount, mmpi_npey)
942 integer :: yourid,jm,jm2,icount,nsize,ierr,jlev,jlev2,jlat,jlat2
943
944 call utl_tmg_start(152,'low-level--gst_barr')
945 if(mmpi_doBarrier) call rpn_comm_barrier('NS',ierr)
946 call utl_tmg_stop(152)
947
948 call utl_tmg_start(154,'low-level--gst_transpose_MtoLAT')
949
950 !$OMP PARALLEL DO PRIVATE(yourid,jlat,jlat2,jlev,jlev2,icount,jm,jm2)
951 do yourid = 0, (mmpi_npey-1)
952 do jlev = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
953 jlev2 = jlev - gst(gstID)%myLevBeg + 1
954 gd_send(:, :, :, jlev2, yourid+1) = 0.0d0
955 do jlat = gst(gstID)%myLatBeg, gst(gstID)%myLatEnd
956 jlat2 = jlat - gst(gstID)%myLatBeg + 1
957 icount = 0
958 do jm = gst(gstID)%allmBeg(yourid+1), gst(gstID)%allmEnd(yourid+1), gst(gstID)%allmSkip(yourid+1)
959 jm2 = 2*jm+1
960 icount = icount + 1
961 gd_send(icount, 1, jlat2, jlev2, yourid+1) = pgd_in(jm2 , jlat, jlev)
962 gd_send(icount, 2, jlat2, jlev2, yourid+1) = pgd_in(jm2+1, jlat, jlev)
963 enddo
964 enddo
965 enddo
966 enddo
967 !$OMP END PARALLEL DO
968
969 nsize = gst(gstID)%maxmCount * 2 * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPEmax
970 if(mmpi_npey.gt.1) then
971 call rpn_comm_alltoall(gd_send, nsize, pre_specTransMpiReal, &
972 gd_recv, nsize, pre_specTransMpiReal, 'NS', ierr)
973 else
974 gd_recv(:,:,:,:,1) = gd_send(:,:,:,:,1)
975 endif
976
977 !$OMP PARALLEL DO PRIVATE(yourid,jlat,jlat2,jlev,jlev2,icount,jm,jm2)
978 do yourid = 0, (mmpi_npey-1)
979 do jlev = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
980 jlev2 = jlev - gst(gstID)%myLevBeg + 1
981 do jlat = gst(gstID)%allLatBeg(yourid+1), gst(gstID)%allLatEnd(yourid+1)
982 jlat2 = jlat - gst(gstID)%allLatBeg(yourid+1) + 1
983 icount = 0
984 do jm = gst(gstID)%mymBeg, gst(gstID)%mymEnd, gst(gstID)%mymSkip
985 jm2 = 2*gst(gstID)%mymIndex(jm)
986 icount = icount + 1
987 pgd_out(jm2-1,jlat,jlev) = gd_recv(icount,1,jlat2,jlev2,yourid+1)
988 pgd_out(jm2 ,jlat,jlev) = gd_recv(icount,2,jlat2,jlev2,yourid+1)
989 enddo
990 enddo
991 enddo
992 enddo
993 !$OMP END PARALLEL DO
994
995 call utl_tmg_stop(154)
996
997 end subroutine transpose2d_LattoM
998
999
1000 subroutine transpose2d_LattoM_kij(pgd_in,pgd_out)
1001 implicit none
1002
1003 ! Arguments:
1004 real(8), intent(in) :: pgd_in(gst(gstID)%maxMyLevCount, gst(gstID)%ni, &
1005 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1006 real(8), intent(out) :: pgd_out(gst(gstID)%maxMyLevCount, 2*gst(gstID)%maxmCount, &
1007 gst(gstID)%nj)
1008
1009 ! Locals:
1010 real(pre_specTransReal), allocatable, save :: gd_send(:,:,:,:,:)
1011 real(pre_specTransReal), allocatable, save :: gd_recv(:,:,:,:,:)
1012 integer :: yourid,jm,jm2,icount,nsize,ierr,jlev,jlat,jlat2
1013
1014 call utl_reAllocate(gd_send, gst(gstID)%maxMyLevCount, gst(gstID)%maxmCount, 2, gst(gstID)%latPerPEmax, mmpi_npey)
1015 call utl_reAllocate(gd_recv, gst(gstID)%maxMyLevCount, gst(gstID)%maxmCount, 2, gst(gstID)%latPerPEmax, mmpi_npey)
1016
1017 call utl_tmg_start(152,'low-level--gst_barr')
1018 if(mmpi_doBarrier) call rpn_comm_barrier('NS',ierr)
1019 call utl_tmg_stop(152)
1020
1021 call utl_tmg_start(154,'low-level--gst_transpose_MtoLAT')
1022
1023 !$OMP PARALLEL DO PRIVATE(yourid,jlat,jlat2,jlev,icount,jm,jm2)
1024 do yourid = 0, (mmpi_npey-1)
1025 gd_send(:, :, :, :, yourid+1) = 0.0d0
1026 do jlat = gst(gstID)%myLatBeg, gst(gstID)%myLatEnd
1027 jlat2 = jlat - gst(gstID)%myLatBeg + 1
1028 icount = 0
1029 do jm = gst(gstID)%allmBeg(yourid+1), gst(gstID)%allmEnd(yourid+1), gst(gstID)%allmSkip(yourid+1)
1030 jm2 = 2*jm+1
1031 icount = icount + 1
1032 do jlev = 1, gst(gstID)%maxMyLevCount
1033 gd_send(jlev, icount, 1, jlat2, yourid+1) = pgd_in(jlev, jm2 , jlat)
1034 gd_send(jlev, icount, 2, jlat2, yourid+1) = pgd_in(jlev, jm2+1, jlat)
1035 enddo
1036 enddo
1037 enddo
1038 enddo
1039 !$OMP END PARALLEL DO
1040
1041 nsize = gst(gstID)%maxmCount * 2 * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPEmax
1042 if(mmpi_npey.gt.1) then
1043 call rpn_comm_alltoall(gd_send, nsize, pre_specTransMpiReal, &
1044 gd_recv, nsize, pre_specTransMpiReal, 'NS', ierr)
1045 else
1046 gd_recv(:,:,:,:,1) = gd_send(:,:,:,:,1)
1047 endif
1048
1049 !$OMP PARALLEL DO PRIVATE(yourid,jlat,jlat2,jlev,icount,jm,jm2)
1050 do yourid = 0, (mmpi_npey-1)
1051 do jlat = gst(gstID)%allLatBeg(yourid+1), gst(gstID)%allLatEnd(yourid+1)
1052 jlat2 = jlat - gst(gstID)%allLatBeg(yourid+1) + 1
1053 icount = 0
1054 do jm = gst(gstID)%mymBeg, gst(gstID)%mymEnd, gst(gstID)%mymSkip
1055 jm2 = 2*gst(gstID)%mymIndex(jm)
1056 icount = icount + 1
1057 do jlev = 1, gst(gstID)%maxMyLevCount
1058 pgd_out(jlev, jm2-1, jlat) = gd_recv(jlev, icount, 1, jlat2, yourid+1)
1059 pgd_out(jlev, jm2 , jlat) = gd_recv(jlev, icount, 2, jlat2, yourid+1)
1060 enddo
1061 enddo
1062 enddo
1063 enddo
1064 !$OMP END PARALLEL DO
1065
1066 call utl_tmg_stop(154)
1067
1068 end subroutine transpose2d_LattoM_kij
1069
1070
1071 subroutine transpose2d_LevtoLon(pgd_in,pgd_out)
1072 implicit none
1073
1074 ! Arguments:
1075 real(8), intent(in) :: pgd_in(gst(gstID)%ni, gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, &
1076 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1077 real(8), intent(out) :: pgd_out(gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1078 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%nk)
1079
1080 ! Locals:
1081 real(pre_specTransReal) :: gd_send(gst(gstID)%lonPerPEmax, gst(gstID)%latPerPEmax, &
1082 gst(gstID)%maxMyLevCount, mmpi_npex)
1083 real(pre_specTransReal) :: gd_recv(gst(gstID)%lonPerPEmax, gst(gstID)%latPerPEmax, &
1084 gst(gstID)%maxMyLevCount, mmpi_npex)
1085 integer :: youridP1,nsize,ierr,jlev,jlev2
1086
1087 call utl_tmg_start(152,'low-level--gst_barr')
1088 if(mmpi_doBarrier) call rpn_comm_barrier('EW',ierr)
1089 call utl_tmg_stop(152)
1090
1091 call utl_tmg_start(155,'low-level--gst_transpose_LEVtoLON')
1092
1093 !$OMP PARALLEL DO PRIVATE(youridP1,jlev,jlev2)
1094 do jlev = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
1095 jlev2 = jlev - gst(gstID)%myLevBeg + 1
1096 do youridP1 = 1, mmpi_npex
1097 gd_send(:, :, jlev2, youridP1) = 0.0d0
1098 gd_send(1:gst(gstID)%allLonPerPE(youridP1), 1:gst(gstID)%latPerPE, jlev2, youridP1) = &
1099 pgd_in(gst(gstID)%allLonBeg(youridP1):gst(gstID)%allLonEnd(youridP1), :, jlev)
1100 enddo
1101 enddo
1102 !$OMP END PARALLEL DO
1103
1104 nsize = gst(gstID)%lonPerPEmax * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPEmax
1105 if(mmpi_npex.gt.1) then
1106 call rpn_comm_alltoall(gd_send,nsize,pre_specTransMpiReal, &
1107 gd_recv,nsize,pre_specTransMpiReal,'EW',ierr)
1108 else
1109 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
1110 endif
1111
1112 !$OMP PARALLEL DO PRIVATE(youridP1,jlev,jlev2)
1113 do youridP1 = 1, mmpi_npex
1114 do jlev=gst(gstID)%allLevBeg(youridP1),gst(gstID)%allLevEnd(youridP1)
1115 jlev2=jlev-gst(gstID)%allLevBeg(youridP1)+1
1116 pgd_out(:, :, jlev) = gd_recv(1:gst(gstID)%lonPerPE, 1:gst(gstID)%latPerPE, jlev2, youridP1)
1117 enddo
1118 enddo
1119 !$OMP END PARALLEL DO
1120
1121 call utl_tmg_stop(155)
1122
1123 end subroutine transpose2d_LevtoLon
1124
1125
1126 subroutine transpose2d_LevtoLon_kij_mpitypes8(pgd_in,pgd_out)
1127 implicit none
1128
1129 ! Arguments:
1130 real(8), intent(in) :: pgd_in(gst(gstID)%maxMyLevCount, gst(gstID)%ni, &
1131 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1132 real(8), intent(out) :: pgd_out(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1133 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1134
1135 ! Locals:
1136 integer :: nsize,ierr
1137
1138 call utl_tmg_start(152,'low-level--gst_barr')
1139 if(mmpi_doBarrier) call rpn_comm_barrier('EW',ierr)
1140 call utl_tmg_stop(152)
1141
1142 call utl_tmg_start(155,'low-level--gst_transpose_LEVtoLON')
1143
1144 nsize = gst(gstID)%lonPerPE * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPE
1145 if(mmpi_npex.gt.1) then
1146 call mpi_alltoall(pgd_in, 1, gst(gstID)%sendType_LevToLon, &
1147 pgd_out, 1, gst(gstID)%recvType_LevToLon, mmpi_comm_EW, ierr)
1148 else
1149 pgd_out(:,:,:) = pgd_in(:,:,:)
1150 endif
1151
1152 call utl_tmg_stop(155)
1153
1154 end subroutine transpose2d_LevtoLon_kij_mpitypes8
1155
1156
1157 subroutine transpose2d_LevtoLon_kij_mpitypes4(pgd_in,pgd_out)
1158 implicit none
1159
1160 ! Arguments:
1161 real(8), intent(in) :: pgd_in(gst(gstID)%maxMyLevCount, gst(gstID)%ni, &
1162 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1163 real(8), intent(out) :: pgd_out(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1164 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1165
1166 ! Locals:
1167 integer :: nsize,ierr
1168 real(4) :: pgd_in_r4(gst(gstID)%maxMyLevCount, gst(gstID)%ni, &
1169 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1170 real(4) :: pgd_out_r4(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1171 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1172
1173 call utl_tmg_start(152,'low-level--gst_barr')
1174 if(mmpi_doBarrier) call rpn_comm_barrier('EW',ierr)
1175 call utl_tmg_stop(152)
1176
1177 call utl_tmg_start(155,'low-level--gst_transpose_LEVtoLON')
1178
1179 nsize = gst(gstID)%lonPerPE * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPE
1180 if(mmpi_npex.gt.1) then
1181 pgd_in_r4(:,:,:) = pgd_in(:,:,:)
1182 call mpi_alltoall(pgd_in_r4, 1, gst(gstID)%sendType_LevToLon, &
1183 pgd_out_r4, 1, gst(gstID)%recvType_LevToLon, mmpi_comm_EW, ierr)
1184 pgd_out(:,:,:) = pgd_out_r4(:,:,:)
1185 else
1186 pgd_out(:,:,:) = pgd_in(:,:,:)
1187 endif
1188
1189 call utl_tmg_stop(155)
1190
1191 end subroutine transpose2d_LevtoLon_kij_mpitypes4
1192
1193
1194 subroutine transpose2d_LevtoLon_kij(pgd_in,pgd_out)
1195 implicit none
1196
1197 ! Arguments:
1198 real(8), intent(in) :: pgd_in(gst(gstID)%maxMyLevCount, gst(gstID)%ni, &
1199 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1200 real(8), intent(out) :: pgd_out(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1201 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1202
1203 ! Locals:
1204 real(pre_specTransReal) :: gd_send(gst(gstID)%maxMyLevCount, gst(gstID)%lonPerPEmax,&
1205 gst(gstID)%latPerPEmax, mmpi_npex)
1206 real(pre_specTransReal) :: gd_recv(gst(gstID)%maxMyLevCount, gst(gstID)%lonPerPEmax,&
1207 gst(gstID)%latPerPEmax, mmpi_npex)
1208 integer :: youridP1, nsize, ierr, yourNumLev
1209
1210 call utl_tmg_start(152,'low-level--gst_barr')
1211 if(mmpi_doBarrier) call rpn_comm_barrier('EW',ierr)
1212 call utl_tmg_stop(152)
1213
1214 call utl_tmg_start(155,'low-level--gst_transpose_LEVtoLON')
1215
1216 !$OMP PARALLEL DO PRIVATE(youridP1)
1217 do youridP1 = 1, mmpi_npex
1218 gd_send(:, :, :, youridP1) = 0.0d0
1219 gd_send(:, 1:gst(gstID)%allLonPerPE(youridP1), 1:gst(gstID)%latPerPE, youridP1) = &
1220 pgd_in(:, gst(gstID)%allLonBeg(youridP1):gst(gstID)%allLonEnd(youridP1), :)
1221 enddo
1222 !$OMP END PARALLEL DO
1223
1224 nsize = gst(gstID)%lonPerPEmax * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPEmax
1225 if(mmpi_npex.gt.1) then
1226 call rpn_comm_alltoall(gd_send, nsize, pre_specTransMpiReal, &
1227 gd_recv, nsize, pre_specTransMpiReal, 'EW', ierr)
1228 else
1229 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
1230 endif
1231
1232 !$OMP PARALLEL DO PRIVATE(youridP1,yourNumLev)
1233 do youridP1 = 1, mmpi_npex
1234 yourNumLev = gst(gstID)%allLevEnd(youridP1) - gst(gstID)%allLevBeg(youridP1) + 1
1235 pgd_out(gst(gstID)%allLevBeg(youridP1):gst(gstID)%allLevEnd(youridP1), :, :) = &
1236 gd_recv(1:yourNumLev, 1:gst(gstID)%lonPerPE, 1:gst(gstID)%latPerPE, youridP1)
1237 enddo
1238 !$OMP END PARALLEL DO
1239
1240 call utl_tmg_stop(155)
1241
1242 end subroutine transpose2d_LevtoLon_kij
1243
1244
1245 subroutine transpose2d_LontoLev(pgd_in,pgd_out)
1246 implicit none
1247
1248 ! Arguments:
1249 real(8), intent(in) :: pgd_in(gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1250 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%nk)
1251 real(8), intent(out) :: pgd_out(gst(gstID)%ni, gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, &
1252 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1253
1254 ! Locals:
1255 real(pre_specTransReal) :: gd_send(gst(gstID)%lonPerPEmax, gst(gstID)%latPerPEmax, &
1256 gst(gstID)%maxMyLevCount, mmpi_npex)
1257 real(pre_specTransReal) :: gd_recv(gst(gstID)%lonPerPEmax, gst(gstID)%latPerPEmax, &
1258 gst(gstID)%maxMyLevCount, mmpi_npex)
1259 integer :: youridP1,nsize,ierr,jlev,jlev2
1260
1261 call utl_tmg_start(152,'low-level--gst_barr')
1262 if(mmpi_doBarrier) call rpn_comm_barrier('EW',ierr)
1263 call utl_tmg_stop(152)
1264
1265 call utl_tmg_start(155,'low-level--gst_transpose_LEVtoLON')
1266
1267 !$OMP PARALLEL DO PRIVATE(youridP1,jlev,jlev2)
1268 do youridP1 = 1, mmpi_npex
1269 do jlev=gst(gstID)%allLevBeg(youridP1),gst(gstID)%allLevEnd(youridP1)
1270 jlev2=jlev-gst(gstID)%allLevBeg(youridP1)+1
1271 gd_send(:, :, jlev2, youridP1) = 0.0d0
1272 gd_send(1:gst(gstID)%lonPerPE, 1:gst(gstID)%latPerPE, jlev2, youridP1) = &
1273 pgd_in(:,:,jlev)
1274 enddo
1275 enddo
1276 !$OMP END PARALLEL DO
1277
1278 nsize = gst(gstID)%lonPerPEmax * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPEmax
1279 if(mmpi_npex.gt.1) then
1280 call rpn_comm_alltoall(gd_send, nsize, pre_specTransMpiReal, &
1281 gd_recv, nsize, pre_specTransMpiReal, 'EW', ierr)
1282 else
1283 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
1284 endif
1285
1286 !$OMP PARALLEL DO PRIVATE(youridP1,jlev,jlev2)
1287 do jlev = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
1288 jlev2 = jlev - gst(gstID)%myLevBeg + 1
1289 do youridP1 = 1, mmpi_npex
1290 pgd_out(gst(gstID)%allLonBeg(youridP1):gst(gstID)%allLonEnd(youridP1), :, jlev) = &
1291 gd_recv(1:gst(gstID)%allLonPerPE(youridP1),1:gst(gstID)%latPerPE,jlev2,youridP1)
1292 enddo
1293 enddo
1294 !$OMP END PARALLEL DO
1295
1296 call utl_tmg_stop(155)
1297
1298 end subroutine transpose2d_LontoLev
1299
1300
1301 subroutine transpose2d_LontoLev_kij_mpitypes8(pgd_in,pgd_out)
1302 implicit none
1303
1304 ! Arguments:
1305 real(8), intent(in) :: pgd_in(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1306 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1307 real(8), intent(out) :: pgd_out(gst(gstID)%maxMyLevCount, gst(gstID)%ni, &
1308 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1309
1310 ! Locals:
1311 integer :: nsize, ierr
1312
1313 call utl_tmg_start(152,'low-level--gst_barr')
1314 if(mmpi_doBarrier) call rpn_comm_barrier('EW',ierr)
1315 call utl_tmg_stop(152)
1316
1317 call utl_tmg_start(155,'low-level--gst_transpose_LEVtoLON')
1318
1319 nsize = gst(gstID)%lonPerPE * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPE
1320 if(mmpi_npex.gt.1) then
1321 call mpi_alltoall(pgd_in, 1, gst(gstID)%sendType_LonToLev, &
1322 pgd_out, 1, gst(gstID)%recvType_LonToLev, mmpi_comm_EW, ierr)
1323 else
1324 pgd_out(:,:,:) = pgd_in(:,:,:)
1325 endif
1326
1327 call utl_tmg_stop(155)
1328
1329 end subroutine transpose2d_LontoLev_kij_mpitypes8
1330
1331
1332 subroutine transpose2d_LontoLev_kij_mpitypes4(pgd_in,pgd_out)
1333 implicit none
1334
1335 ! Arguments:
1336 real(8), intent(in) :: pgd_in(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1337 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1338 real(8), intent(out) :: pgd_out(gst(gstID)%maxMyLevCount, gst(gstID)%ni, &
1339 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1340
1341 ! Locals:
1342 integer :: nsize, ierr
1343 real(4) :: pgd_in_r4(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1344 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1345 real(4) :: pgd_out_r4(gst(gstID)%maxMyLevCount, gst(gstID)%ni, &
1346 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1347
1348 call utl_tmg_start(152,'low-level--gst_barr')
1349 if(mmpi_doBarrier) call rpn_comm_barrier('EW',ierr)
1350 call utl_tmg_stop(152)
1351
1352 call utl_tmg_start(155,'low-level--gst_transpose_LEVtoLON')
1353
1354 nsize = gst(gstID)%lonPerPE * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPE
1355 if(mmpi_npex.gt.1) then
1356 pgd_in_r4(:,:,:) = pgd_in(:,:,:)
1357 call mpi_alltoall(pgd_in_r4, 1, gst(gstID)%sendType_LonToLev, &
1358 pgd_out_r4, 1, gst(gstID)%recvType_LonToLev, mmpi_comm_EW, ierr)
1359 pgd_out(:,:,:) = pgd_out_r4(:,:,:)
1360 else
1361 pgd_out(:,:,:) = pgd_in(:,:,:)
1362 endif
1363
1364 call utl_tmg_stop(155)
1365
1366 end subroutine transpose2d_LontoLev_kij_mpitypes4
1367
1368
1369 subroutine transpose2d_LontoLev_kij(pgd_in,pgd_out)
1370 implicit none
1371
1372 ! Arguments:
1373 real(8), intent(in) :: pgd_in(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1374 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1375 real(8), intent(out) :: pgd_out(gst(gstID)%maxMyLevCount, gst(gstID)%ni, &
1376 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
1377
1378 ! Locals:
1379 real(pre_specTransReal) :: gd_send(gst(gstID)%maxMyLevCount, gst(gstID)%lonPerPEmax,&
1380 gst(gstID)%latPerPEmax, mmpi_npex)
1381 real(pre_specTransReal) :: gd_recv(gst(gstID)%maxMyLevCount, gst(gstID)%lonPerPEmax,&
1382 gst(gstID)%latPerPEmax, mmpi_npex)
1383 integer :: youridP1,nsize,ierr,jlev,jlev2,yourNumLev
1384
1385 call utl_tmg_start(152,'low-level--gst_barr')
1386 if(mmpi_doBarrier) call rpn_comm_barrier('EW',ierr)
1387 call utl_tmg_stop(152)
1388
1389 call utl_tmg_start(155,'low-level--gst_transpose_LEVtoLON')
1390
1391 !$OMP PARALLEL DO PRIVATE(youridP1,yourNumLev)
1392 do youridP1 = 1, mmpi_npex
1393 yourNumLev = gst(gstID)%allLevEnd(youridP1) - gst(gstID)%allLevBeg(youridP1) + 1
1394 gd_send(:, :, :, youridP1) = 0.0d0
1395 gd_send(1:yourNumLev, 1:gst(gstID)%lonPerPE, 1:gst(gstID)%latPerPE, youridP1) = &
1396 pgd_in(gst(gstID)%allLevBeg(youridP1):gst(gstID)%allLevEnd(youridP1),:,:)
1397 enddo
1398 !$OMP END PARALLEL DO
1399
1400 nsize = gst(gstID)%lonPerPEmax * gst(gstID)%maxMyLevCount * gst(gstID)%latPerPEmax
1401 if(mmpi_npex.gt.1) then
1402 call rpn_comm_alltoall(gd_send, nsize, pre_specTransMpiReal, &
1403 gd_recv, nsize, pre_specTransMpiReal, 'EW', ierr)
1404 else
1405 gd_recv(:,:,:,1) = gd_send(:,:,:,1)
1406 endif
1407
1408 !$OMP PARALLEL DO PRIVATE(youridP1,jlev,jlev2)
1409 do youridP1 = 1, mmpi_npex
1410 pgd_out(:, gst(gstID)%allLonBeg(youridP1):gst(gstID)%allLonEnd(youridP1), :) = &
1411 gd_recv(:, 1:gst(gstID)%allLonPerPE(youridP1), 1:gst(gstID)%latPerPE, youridP1)
1412 enddo
1413 !$OMP END PARALLEL DO
1414
1415 call utl_tmg_stop(155)
1416
1417 end subroutine transpose2d_LontoLev_kij
1418
1419!-------------------------------------------------------------------------------
1420! Subroutines to re-order the u and v wind components for mpi version of spgd
1421! and spgda
1422!-------------------------------------------------------------------------------
1423
1424 subroutine interleaveWinds_sp(psp,nflev)
1425 implicit none
1426
1427 ! Arguments:
1428 integer, intent(in) :: nflev
1429 real(8), intent(inout) :: psp(gst(gstID)%myNla,2,gst(gstID)%nk)
1430
1431 ! Locals:
1432 real(8) :: tempvalues(2,nflev*2)
1433 integer :: jk, ila
1434
1435 !$OMP PARALLEL DO PRIVATE (ILA,JK,TEMPVALUES)
1436 do ila = 1, gst(gstID)%myNla
1437 do jk = 1, nflev
1438 ! place u in new position in temporary array
1439 tempvalues(:,(jk*2)-1) = psp(ila,:,jk)
1440 ! place v in new position in temporary array
1441 tempvalues(:,jk*2) = psp(ila,:,jk+nflev)
1442 enddo
1443 ! move contents of temporary array back to original array
1444 psp(ila,:,1:2*nflev) = tempvalues(:,1:2*nflev)
1445 enddo
1446 !$OMP END PARALLEL DO
1447
1448 end subroutine interleaveWinds_sp
1449
1450
1451 subroutine unInterleaveWinds_sp(psp,nflev)
1452 implicit none
1453
1454 ! Arguments:
1455 integer, intent(in) :: nflev
1456 real(8), intent(inout) :: psp(gst(gstID)%myNla,2,gst(gstID)%nk)
1457
1458 ! Locals:
1459 real(8) :: tempvalues(2,nflev*2)
1460 integer :: jk, ila
1461
1462 !$OMP PARALLEL DO PRIVATE (ILA,JK,TEMPVALUES)
1463 do ila = 1, gst(gstID)%myNla
1464 do jk = 1, nflev
1465 ! place u in original position in temporary array
1466 tempvalues(:,jk) = psp(ila,:,(jk*2)-1)
1467 ! place v in original position in temporary array
1468 tempvalues(:,jk+nflev) = psp(ila,:,jk*2)
1469 enddo
1470 ! move contents of temporary array back to original array
1471 psp(ila,:,1:2*nflev) = tempvalues(:,1:2*nflev)
1472 enddo
1473 !$OMP END PARALLEL DO
1474
1475 end subroutine unInterleaveWinds_sp
1476
1477
1478 subroutine interleaveWinds_gd(pgd,nflev)
1479 implicit none
1480
1481 ! Arguments:
1482 integer, intent(in) :: nflev
1483 real(8), intent(inout) :: pgd(gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1484 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%nk)
1485
1486 ! Locals:
1487 real(8) :: tempvalues(nflev*2)
1488 integer :: jlat, jk, jlon
1489
1490 !$OMP PARALLEL DO PRIVATE (JLAT,JLON,JK,TEMPVALUES)
1491 do jlat = gst(gstID)%myLatBeg, gst(gstID)%myLatEnd
1492 do jlon = gst(gstID)%myLonBeg, gst(gstID)%myLonEnd
1493 do jk = 1, nflev
1494 ! place u in original position in temporary array
1495 tempvalues((jk*2)-1) = pgd(jlon,jlat,jk)
1496 ! place v in original position in temporary array
1497 tempvalues(jk*2) = pgd(jlon,jlat,jk+nflev)
1498 enddo
1499 ! move contents of temporary array back to original array
1500 pgd(jlon,jlat,1:2*nflev) = tempvalues(1:2*nflev)
1501 enddo
1502 enddo
1503 !$OMP END PARALLEL DO
1504
1505 end subroutine interleaveWinds_gd
1506
1507
1508 subroutine unInterleaveWinds_gd(pgd,nflev)
1509 implicit none
1510
1511 ! Arguments:
1512 integer, intent(in) :: nflev
1513 real(8), intent(inout) :: pgd(gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1514 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%nk)
1515
1516 ! Locals:
1517 real(8) :: tempvalues(nflev*2)
1518 integer :: jlat, jk, jlon
1519
1520 !$OMP PARALLEL DO PRIVATE (JLAT,JLON,JK,TEMPVALUES)
1521 do jlat = gst(gstID)%myLatBeg, gst(gstID)%myLatEnd
1522 do jlon = gst(gstID)%myLonBeg, gst(gstID)%myLonEnd
1523 do jk = 1, nflev
1524 ! place u in original position in temporary array
1525 tempvalues(jk) = pgd(jlon,jlat,(jk*2)-1)
1526 ! place v in original position in temporary array
1527 tempvalues(jk+nflev) = pgd(jlon,jlat,jk*2)
1528 enddo
1529 ! move contents of temporary array back to original array
1530 pgd(jlon,jlat,1:2*nflev) = tempvalues(1:2*nflev)
1531 enddo
1532 enddo
1533 !$OMP END PARALLEL DO
1534
1535 end subroutine unInterleaveWinds_gd
1536
1537!-------------------------------------------------------------------------------
1538! Main spectral transform subroutines
1539!-------------------------------------------------------------------------------
1540
1541 subroutine gst_spgd(psp,pgd,nflev)
1542 implicit none
1543
1544 ! Arguments:
1545 integer, intent(in) :: nflev
1546 real(8), intent(inout) :: psp(gst(gstID)%myNla,2,gst(gstID)%nk)
1547 real(8), intent(out) :: pgd(gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1548 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%nk)
1549
1550 ! Locals:
1551 real(8), allocatable :: psp2(:,:,:),pgd2(:,:,:),pgd3(:,:,:)
1552 integer :: jlat, jk, jlon
1553
1554 ! check if this mpi task will deal with winds during Legendre transform
1555 if(gst(gstID)%myLevBeg.le.2*nflev) then
1556 ! ensure that the number of levels on this mpi task is even to allow interleaving of u and v
1557 ! only necessary when number of levels on an mpi task is less than all wind levels (2nd condition)
1558 if( (mod(gst(gstID)%myLevCount,2).ne.0) .and. (gst(gstID)%myLevCount < 2*nflev) ) then
1559 write(*,*) 'GST_SPGD: myLevCount = ',gst(gstID)%myLevCount
1560 call utl_abort('GST_SPGD: Number of levels on this mpi task must be even!')
1561 endif
1562 endif
1563
1564 allocate(psp2(gst(gstID)%nla, 2, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
1565 allocate(pgd2(2*gst(gstID)%maxmcount, gst(gstID)%nj, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
1566 allocate(pgd3(gst(gstID)%ni, gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
1567
1568 ! 1.0 First reorder wind components to have u and v for same level on same mpi task
1569 call interleaveWinds_sp(psp,nflev)
1570
1571 ! 1.1 Transpose data along npex from N to Levels
1572 call transpose2d_NtoLev(psp,psp2)
1573
1574 ! 1.2 Inverse Legendre transform
1575 call utl_tmg_start(150,'low-level--gst_lt')
1576 call spgdpar(psp2,pgd2,nflev)
1577 call utl_tmg_stop(150)
1578 deallocate(psp2)
1579
1580 ! 1.3 Transpose data along npey from M to Latitudes
1581 call transpose2d_MtoLat(pgd2,pgd3)
1582 deallocate(pgd2)
1583
1584 !$OMP PARALLEL DO PRIVATE (JLAT,JLON,JK)
1585 ! 2.1 Reset to zero the modes that are not part of the truncation
1586 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
1587 do jlat = gst(gstID)%myLatBeg, gst(gstID)%myLatEnd
1588 do jlon = 2*(gst(gstID)%ntrunc+1)+1, gst(gstID)%ni
1589 pgd3(jlon,jlat,jk) = 0.d0
1590 enddo
1591 enddo
1592 enddo
1593 !$OMP END PARALLEL DO
1594
1595 ! 2.2 Apply the FFT
1596 call utl_tmg_start(151,'low-level--gst_fft')
1597 call fft3dvar(pgd3,+1)
1598 call utl_tmg_stop(151)
1599
1600 ! 2.3 Transpose data along npex from Levels to Longitudes
1601 call transpose2d_LevtoLon(pgd3,pgd)
1602 deallocate(pgd3)
1603
1604 ! 2.4 Now undo reordering of wind components
1605 call unInterleaveWinds_gd(pgd,nflev)
1606
1607 end subroutine gst_spgd
1608
1609
1610 subroutine gst_gdsp(psp,pgd,nflev)
1611 implicit none
1612
1613 ! Arguments:
1614 integer, intent(in) :: nflev
1615 real(8), intent(out) :: psp(gst(gstID)%myNla,2,gst(gstID)%nk)
1616 real(8), intent(inout) :: pgd(gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1617 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd,gst(gstID)%nk)
1618
1619 ! Locals:
1620 real(8), allocatable :: psp2(:,:,:),pgd2(:,:,:),pgd3(:,:,:)
1621
1622 ! check if this mpi task will deal with winds during Legendre transform
1623 if(gst(gstID)%myLevBeg.le.2*nflev) then
1624 ! ensure that the number of levels on this mpi task is even to allow interleaving of u and v
1625 ! only necessary when number of levels on an mpi task is less than all wind levels (2nd condition)
1626 if( (mod(gst(gstID)%myLevCount,2).ne.0) .and. (gst(gstID)%myLevCount < 2*nflev) ) then
1627 write(*,*) 'GST_GDSP: myLevCount = ',gst(gstID)%myLevCount
1628 call utl_abort('GST_GDSP: Number of levels on this mpi task must be even!')
1629 endif
1630 endif
1631
1632 allocate(psp2(gst(gstID)%nla, 2, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
1633 allocate(pgd2(2*gst(gstID)%maxmcount, gst(gstID)%nj, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
1634 allocate(pgd3(gst(gstID)%ni, gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
1635
1636 ! 1.0 First reorder wind components to have u and v for same level on same mpi task
1637 call interleaveWinds_gd(pgd,nflev)
1638
1639 ! Transpose data along npex from Longitudes to Levels
1640 call transpose2d_LontoLev(pgd,pgd3)
1641
1642 ! 1. Fourier transform all fields for all latitudes
1643 call utl_tmg_start(151,'low-level--gst_fft')
1644 call fft3dvar(pgd3,-1)
1645 call utl_tmg_stop(151)
1646
1647 ! Transpose data along npey from Latitudes to M
1648 call transpose2d_LattoM(pgd3,pgd2)
1649 deallocate(pgd3)
1650
1651 ! 2. Direct Legendre transform including wind transformations
1652 call utl_tmg_start(150,'low-level--gst_lt')
1653 call gdsppar(psp2,pgd2,nflev)
1654 call utl_tmg_stop(150)
1655 deallocate(pgd2)
1656
1657 ! Transpose data along npex from Levels to N
1658 call transpose2d_LevtoN(psp2,psp)
1659 deallocate(psp2)
1660
1661 ! 2.4 Now undo reordering of wind components
1662 call unInterleaveWinds_sp(psp,nflev)
1663
1664 end subroutine gst_gdsp
1665
1666
1667 subroutine spgdpar(psp,pgd2,nflev)
1668 !
1669 !:Purpose: Inverse spectral transform(PARALLEL LOOP)
1670 !
1671 implicit none
1672
1673 ! Arguments:
1674 integer, intent(in) :: nflev
1675 real(8), intent(in) :: psp(gst(gstID)%nla,2,gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1676 real(8), intent(out) :: pgd2(2*gst(gstID)%maxmcount,gst(gstID)%nj,&
1677 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1678
1679 ! Locals
1680 integer :: jj, jj2, jm, jn, ilonr, iloni, jk, jk2, ila, inm
1681 real(8) :: zjm, factor
1682 real(8) :: dlalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
1683 real(8) :: dldalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
1684 real(8) :: zfms(gst(gstID)%njlath+1,2,gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1685 real(8) :: zfma(gst(gstID)%njlath+1,2,gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1686 real(8) :: dlsp(0:gst(gstID)%ntrunc,2,gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1687
1688 ! Inverse Legendre transform
1689
1690 !$OMP PARALLEL DO PRIVATE(DLALP,DLDALP,DLSP,ZFMS,ZFMA, &
1691 !$OMP INM,ILA,JM,JN,JK,JK2,JJ,JJ2,ZJM,ILONR,ILONI,FACTOR)
1692 do jm = gst(gstID)%mymBeg, gst(gstID)%mymEnd, gst(gstID)%mymSkip
1693
1694 ilonr = 2*gst(gstID)%mymIndex(jm)-1
1695 iloni = 2*gst(gstID)%mymIndex(jm)
1696 zjm = real(jm,8)
1697
1698 ! 2.1 Copy global spectral state into local spectral state
1699 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
1700 do jn = jm, gst(gstID)%ntrunc
1701 ila = gst(gstID)%nind(jm) + jn - jm
1702 inm = jn - jm
1703 if(jk > 2*nflev) then
1704 ! Scalar fields
1705 dlsp(inm,1,jk) = psp(ila,1,jk)
1706 dlsp(inm,2,jk) = psp(ila,2,jk)
1707 else
1708 ! Vector fields
1709 dlsp(inm,1,jk) = psp(ila,1,jk)*gst(gstID)%r1snp1(ila)
1710 dlsp(inm,2,jk) = psp(ila,2,jk)*gst(gstID)%r1snp1(ila)
1711 endif
1712 enddo
1713 enddo
1714
1715 ! 2.2 Get Legendre polynomial (and its derivative) for all latitudes
1716 ! but for the chosen value of "m" from the global array
1717 call getalp (dlalp,dldalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc,jm)
1718
1719 ! 2.3 Perform the inverse Legendre transform for all fields
1720 call leginv2d(jm,zfms,zfma,dlsp,dlalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc)
1721
1722 ! 2.4 Passage to Fourier space
1723
1724 ! Scalar fields
1725 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
1726 if(jk > 2*nflev) then
1727 do jj = 1, gst(gstID)%nj
1728 jj2 = gst(gstID)%nj - jj + 1
1729 if(jj.le.gst(gstID)%njlath) then
1730 pgd2(ilonr,jj2,jk) = zfms(jj,1,jk) + zfma(jj,1,jk)
1731 pgd2(iloni,jj2,jk) = zfms(jj,2,jk) + zfma(jj,2,jk)
1732 else
1733 pgd2(ilonr,jj2,jk) = zfms(jj2,1,jk) - zfma(jj2,1,jk)
1734 pgd2(iloni,jj2,jk) = zfms(jj2,2,jk) - zfma(jj2,2,jk)
1735 endif
1736 enddo
1737 endif
1738 enddo
1739
1740 ! Vector fields: Note that u and v are interleaved in mode 5!
1741 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd, 2
1742 if(jk <= 2*nflev) then
1743 jk2 = jk + 1 ! jk is u, jk2 is v
1744 do jj = 1, gst(gstID)%nj
1745 jj2 = gst(gstID)%nj - jj + 1
1746 if(jj.le.gst(gstID)%njlath) then
1747 pgd2(ilonr,jj2,jk ) = -zjm*(zfms(jj,2,jk2) + zfma(jj,2,jk2))
1748 pgd2(iloni,jj2,jk ) = zjm*(zfms(jj,1,jk2) + zfma(jj,1,jk2))
1749 pgd2(ilonr,jj2,jk2) = -zjm*(zfms(jj,2,jk) + zfma(jj,2,jk))
1750 pgd2(iloni,jj2,jk2) = zjm*(zfms(jj,1,jk) + zfma(jj,1,jk))
1751 else
1752 pgd2(ilonr,jj2,jk) = -zjm*(zfms(jj2,2,jk2) - zfma(jj2,2,jk2))
1753 pgd2(iloni,jj2,jk) = zjm*(zfms(jj2,1,jk2) - zfma(jj2,1,jk2))
1754 pgd2(ilonr,jj2,jk2) = -zjm*(zfms(jj2,2,jk) - zfma(jj2,2,jk))
1755 pgd2(iloni,jj2,jk2) = zjm*(zfms(jj2,1,jk) - zfma(jj2,1,jk))
1756 endif
1757 enddo
1758 endif
1759 enddo
1760
1761 ! 2.5 Completion of the computation of the winds in Fourier space
1762 call leginv2d(jm,zfma,zfms,dlsp,dldalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc)
1763
1764 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd, 2
1765 if(jk <= 2*nflev) then ! only for winds
1766 jk2 = jk + 1
1767 do jj = 1, gst(gstID)%nj
1768 jj2 = gst(gstID)%nj - jj + 1
1769 factor = ec_ra / cos(gst(gstID)%rlati(jj))
1770 if(jj.ne.jj2) then
1771 ! For latitudes not exactly at equator
1772 if(jj.le.gst(gstID)%njlath) then
1773 ! northern latitudes
1774 ! u component
1775 pgd2(ilonr,jj2,jk) = factor * ( pgd2(ilonr,jj2,jk) - (zfms(jj,1,jk) + zfma(jj,1,jk)) )
1776 pgd2(iloni,jj2,jk) = factor * ( pgd2(iloni,jj2,jk) - (zfms(jj,2,jk) + zfma(jj,2,jk)) )
1777 ! v component
1778 pgd2(ilonr,jj2,jk2) = factor * ( pgd2(ilonr,jj2,jk2) + (zfms(jj,1,jk2) + zfma(jj,1,jk2)) )
1779 pgd2(iloni,jj2,jk2) = factor * ( pgd2(iloni,jj2,jk2) + (zfms(jj,2,jk2) + zfma(jj,2,jk2)) )
1780 else
1781 ! southern latitudes
1782 ! u component
1783 pgd2(ilonr,jj2,jk ) = factor * ( pgd2(ilonr,jj2,jk ) -(zfms(jj2,1,jk ) - zfma(jj2,1,jk )) )
1784 pgd2(iloni,jj2,jk ) = factor * ( pgd2(iloni,jj2,jk ) -(zfms(jj2,2,jk ) - zfma(jj2,2,jk )) )
1785 ! v component
1786 pgd2(ilonr,jj2,jk2) = factor * ( pgd2(ilonr,jj2,jk2) +(zfms(jj2,1,jk2) - zfma(jj2,1,jk2)) )
1787 pgd2(iloni,jj2,jk2) = factor * ( pgd2(iloni,jj2,jk2) +(zfms(jj2,2,jk2) - zfma(jj2,2,jk2)) )
1788 endif
1789 else
1790 ! Special case for the equator (jj.eq.jj2)
1791 write(*,*) 'SPGDPAR: special case of jj.eq.jj2!!!'
1792 ! u component
1793 pgd2(ilonr,jj2,jk ) = factor * ( pgd2(ilonr,jj2,jk ) -(zfms(jj,1,jk ) + zfma(jj,1,jk )) )
1794 pgd2(iloni,jj2,jk ) = factor * ( pgd2(iloni,jj2,jk ) -(zfms(jj,2,jk ) + zfma(jj,2,jk )) )
1795 ! v component
1796 pgd2(ilonr,jj2,jk2) = factor * ( pgd2(ilonr,jj2,jk2) +(zfms(jj,1,jk2) + zfma(jj,1,jk2)) )
1797 pgd2(iloni,jj2,jk2) = factor * ( pgd2(iloni,jj2,jk2) +(zfms(jj,2,jk2) + zfma(jj,2,jk2)) )
1798 endif
1799 enddo ! jj
1800 endif
1801 enddo ! jk
1802 enddo ! end loop on m
1803 !$OMP END PARALLEL DO
1804
1805 end subroutine spgdpar
1806
1807
1808 subroutine gdsppar(psp,pgd2,nflev)
1809 implicit none
1810
1811 ! Arguments:
1812 integer, intent(in) :: nflev
1813 real(8), intent(out) :: psp(gst(gstID)%nla,2,gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1814 real(8), intent(in) :: pgd2(2*gst(gstID)%maxmcount,gst(gstID)%nj, &
1815 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1816
1817 ! Locals:
1818 integer :: jj, jj2, jk, jk2, ilonr, iloni, jm ,ila, inm, jn
1819 real(8) :: zjm, dlrwt(gst(gstID)%nj), dlrwocs(gst(gstID)%nj)
1820 real(8) :: dlalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
1821 real(8) :: dldalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
1822 real(8) :: dlsp(0:gst(gstID)%ntrunc,2, &
1823 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1824 real(8) :: dlsp2(0:gst(gstID)%ntrunc,2, &
1825 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1826 real(8) :: zfms(gst(gstID)%njlath+1,2, &
1827 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1828 real(8) :: zfma(gst(gstID)%njlath+1,2, &
1829 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
1830
1831 ! 1. Adjustment needed when an odd number of latitudes is considered
1832 dlrwt(:) = gst(gstID)%rwt(:)
1833 !dlrwocs(:) = gst(gstID)%rwocs(:)
1834 do jj = 1, gst(gstID)%nj
1835 dlrwocs(jj) = gst(gstID)%rwt(jj)/(ec_ra*cos(gst(gstID)%rlati(jj)))
1836 enddo
1837 if (mod(gst(gstID)%nj,2).ne.0) then
1838 dlrwt(gst(gstID)%njlath) = dlrwt(gst(gstID)%njlath)/2.d0
1839 dlrwocs(gst(gstID)%njlath) = dlrwocs(gst(gstID)%njlath)/2.d0
1840 end if
1841
1842 !$OMP PARALLEL DO PRIVATE(DLALP,DLDALP,DLSP,DLSP2,ZFMS,ZFMA, &
1843 !$OMP INM,ILA,JM,JN,JK,JK2,JJ,JJ2,ZJM,ILONR,ILONI)
1844 do jm = gst(gstID)%mymBeg, gst(gstID)%mymEnd, gst(gstID)%mymSkip
1845
1846 ilonr = 2*gst(gstID)%mymIndex(jm)-1
1847 iloni = 2*gst(gstID)%mymIndex(jm)
1848 zjm = real(jm,8)
1849
1850 ! 2.1 Fetch the Legendre functions and their derivatives for this choice of "m"
1851 call getalp(dlalp,dldalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc,jm)
1852
1853 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
1854 do jj = 1, gst(gstID)%njlath
1855 zfms(jj,1,jk) = 0.0d0
1856 zfms(jj,2,jk) = 0.0d0
1857 zfma(jj,1,jk) = 0.0d0
1858 zfma(jj,2,jk) = 0.0d0
1859 enddo
1860 enddo
1861
1862 ! 2.2 Build the symmetric and anti-symmetric Fourier coefficients including
1863 ! the appropriate quadrature weights (see scientific notes)
1864 do jj = 1, gst(gstID)%njlath
1865 jj2 = gst(gstID)%nj - jj + 1
1866
1867 ! 2.2.1 Coefficients for scalar fields
1868 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
1869 if(jk .gt. 2*nflev) then
1870 ! symmetric coefficients
1871 zfms(jj,1,jk) = dlrwt(jj)*(pgd2(ilonr,jj2,jk) + pgd2(ilonr,jj,jk))
1872 zfms(jj,2,jk) = dlrwt(jj)*(pgd2(iloni,jj2,jk) + pgd2(iloni,jj,jk))
1873 ! antisymmetric coefficients
1874 zfma(jj,1,jk) = dlrwt(jj)*(pgd2(ilonr,jj2,jk) - pgd2(ilonr,jj,jk))
1875 zfma(jj,2,jk) = dlrwt(jj)*(pgd2(iloni,jj2,jk) - pgd2(iloni,jj,jk))
1876 endif
1877 enddo
1878
1879 ! 2.2.2 Coefficients associated with the wind fields
1880 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd, 2
1881 if(jk .le. 2*nflev) then
1882 jk2 = jk + 1 ! jk is u, jk2 is v
1883 ! vorticity: symmetric coefficients
1884 zfms(jj,1,jk) = -zjm*dlrwocs(jj)*(pgd2(iloni,jj2,jk2) + pgd2(iloni,jj,jk2))
1885 zfms(jj,2,jk) = zjm*dlrwocs(jj)*(pgd2(ilonr,jj2,jk2) + pgd2(ilonr,jj,jk2))
1886 ! vorticity: antisymmetric coefficients
1887 zfma(jj,1,jk) = -zjm*dlrwocs(jj)*(pgd2(iloni,jj2,jk2) - pgd2(iloni,jj,jk2))
1888 zfma(jj,2,jk) = zjm*dlrwocs(jj)*(pgd2(ilonr,jj2,jk2) - pgd2(ilonr,jj,jk2))
1889 ! divergence: symmetric coefficients
1890 zfms(jj,1,jk2) = -zjm*dlrwocs(jj)*(pgd2(iloni,jj2,jk) + pgd2(iloni,jj,jk))
1891 zfms(jj,2,jk2) = zjm*dlrwocs(jj)*(pgd2(ilonr,jj2,jk) + pgd2(ilonr,jj,jk))
1892 ! divergence: antisymmetric coefficients
1893 zfma(jj,1,jk2) = -zjm*dlrwocs(jj)*(pgd2(iloni,jj2,jk) - pgd2(iloni,jj,jk))
1894 zfma(jj,2,jk2) = zjm*dlrwocs(jj)*(pgd2(ilonr,jj2,jk) - pgd2(ilonr,jj,jk))
1895 endif
1896 enddo
1897 enddo
1898
1899 ! 2.3 First one with ALP for all scalar fields and for half the terms
1900 ! required to define the divergence and vorticity
1901 call legdir2d(jm,zfms,zfma,dlsp,dlalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc)
1902
1903 ! 2.4 Second transform with DALP to complete the construction of the
1904 ! vorticity and divergence fields
1905 do jj = 1, gst(gstID)%njlath
1906 jj2 = gst(gstID)%nj - jj + 1
1907 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd, 2
1908 if(jk <= 2*nflev) then ! only for winds
1909 jk2 = jk + 1
1910
1911 ! symmetric coefficients for zonal wind
1912 zfms(jj,1,jk) = dlrwocs(jj)*(pgd2(ilonr,jj2,jk) + pgd2(ilonr,jj,jk))
1913 zfms(jj,2,jk) = dlrwocs(jj)*(pgd2(iloni,jj2,jk) + pgd2(iloni,jj,jk))
1914 ! antisymmetric coefficients for zonal wind
1915 zfma(jj,1,jk) = dlrwocs(jj)*(pgd2(ilonr,jj2,jk) - pgd2(ilonr,jj,jk))
1916 zfma(jj,2,jk) = dlrwocs(jj)*(pgd2(iloni,jj2,jk) - pgd2(iloni,jj,jk))
1917
1918 ! symmetric coefficients for zonal wind
1919 zfms(jj,1,jk2) = -dlrwocs(jj)*(pgd2(ilonr,jj2,jk2) + pgd2(ilonr,jj,jk2))
1920 zfms(jj,2,jk2) = -dlrwocs(jj)*(pgd2(iloni,jj2,jk2) + pgd2(iloni,jj,jk2))
1921 ! antisymmetric coefficients for zonal wind
1922 zfma(jj,1,jk2) = -dlrwocs(jj)*(pgd2(ilonr,jj2,jk2) - pgd2(ilonr,jj,jk2))
1923 zfma(jj,2,jk2) = -dlrwocs(jj)*(pgd2(iloni,jj2,jk2) - pgd2(iloni,jj,jk2))
1924
1925 endif
1926 enddo
1927 enddo
1928
1929 call legdir2d(jm,zfma,zfms,dlsp2,dldalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc)
1930
1931 ! 2.5 Transfer the result in the global state
1932 do jn = jm, gst(gstID)%ntrunc
1933 ila = gst(gstid)%nind(jm) + jn - jm
1934 inm = jn - jm
1935 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
1936 if(jk <= 2*nflev) then ! for winds
1937 psp(ila,1,jk) = dlsp(inm,1,jk) + dlsp2(inm,1,jk)
1938 psp(ila,2,jk) = dlsp(inm,2,jk) + dlsp2(inm,2,jk)
1939 else ! for scalar fields
1940 psp(ila,1,jk) = dlsp(inm,1,jk)
1941 psp(ila,2,jk) = dlsp(inm,2,jk)
1942 endif
1943 enddo
1944 enddo
1945
1946 enddo ! jm
1947 !$OMP END PARALLEL DO
1948 end subroutine gdsppar
1949
1950
1951 subroutine gst_spgda(psp,pgd,nflev)
1952 implicit none
1953
1954 ! Arguments:
1955 integer, intent(in) :: nflev
1956 real(8), intent(inout) :: psp(gst(gstID)%myNla,2,gst(gstID)%nk)
1957 real(8), intent(inout) :: pgd(gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
1958 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%nk)
1959
1960 ! Locals:
1961 real(8), allocatable :: psp2(:,:,:),pgd2(:,:,:),pgd3(:,:,:)
1962
1963 ! check if this mpi task will deal with winds during Legendre transform
1964 if(gst(gstID)%myLevBeg.le.2*nflev) then
1965 ! ensure that the number of levels on this mpi task is even to allow interleaving of u and v
1966 ! only necessary when number of levels on an mpi task is less than all wind levels (2nd condition)
1967 if( (mod(gst(gstID)%myLevCount,2).ne.0) .and. (gst(gstID)%myLevCount < 2*nflev) ) then
1968 write(*,*) 'GST_SPGDA: myLevCount = ',gst(gstID)%myLevCount
1969 call utl_abort('GST_SPGDA: Number of levels on this mpi task must be even!')
1970 endif
1971 endif
1972
1973 allocate(psp2(gst(gstID)%nla, 2, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
1974 allocate(pgd2(2*gst(gstID)%maxmcount, gst(gstID)%nj, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
1975 allocate(pgd3(gst(gstID)%ni, gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
1976
1977 call adjnorm(pgd)
1978
1979 ! First reorder wind components to have u and v for same level on same mpi task
1980 call interleaveWinds_gd(pgd,nflev)
1981
1982 ! Transpose data along npex from Longitudes to Levels
1983 call transpose2d_LontoLev(pgd,pgd3)
1984
1985 ! Fourier transform all fields for all latitudes
1986 call utl_tmg_start(151,'low-level--gst_fft')
1987 call fft3dvar(pgd3,-1)
1988 call utl_tmg_stop(151)
1989
1990 ! Transpose data along npey from Latitudes to M
1991 call transpose2d_LattoM(pgd3,pgd2)
1992 deallocate(pgd3)
1993
1994 ! Direct Legendre transform including wind transformations
1995 call utl_tmg_start(150,'low-level--gst_lt')
1996 call spgdapar(psp2,pgd2,nflev)
1997 call utl_tmg_stop(150)
1998 deallocate(pgd2)
1999
2000 ! Transpose data along npex from Levels to N
2001 call transpose2d_LevtoN(psp2,psp)
2002 deallocate(psp2)
2003
2004 ! Now undo reordering of wind components
2005 call unInterleaveWinds_sp(psp,nflev)
2006
2007 end subroutine gst_spgda
2008
2009
2010 subroutine spgdapar(psp,pgd2,nflev)
2011 implicit none
2012
2013 ! Arguments:
2014 integer, intent(in) :: nflev
2015 real(8), intent(out) :: psp(gst(gstID)%nla,2,gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2016 real(8), intent(in) :: pgd2(2*gst(gstID)%maxmCount,gst(gstID)%nj, &
2017 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2018
2019 ! Locals:
2020 integer :: jj, jj2, jk, jk2, ilonr, iloni, jm ,ila, inm, jn
2021 real(8) :: zjm,factor,dlrwt(gst(gstID)%nj)
2022 real(8) :: dlalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
2023 real(8) :: dldalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
2024 real(8) :: dlsp(0:gst(gstID)%ntrunc,2, &
2025 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2026 real(8) :: dlsp2(0:gst(gstID)%ntrunc,2, &
2027 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2028 real(8) :: zfms(gst(gstID)%njlath+1,2, &
2029 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2030 real(8) :: zfma(gst(gstID)%njlath+1,2, &
2031 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2032
2033 ! 1. Set up according to the desired grid
2034 ! ---------------------
2035 dlrwt(:) = gst(gstID)%rwt(:)
2036 if (mod(gst(gstID)%nj,2).ne.0) then
2037 dlrwt(gst(gstID)%njlath) = dlrwt(gst(gstID)%njlath)/2.d0
2038 end if
2039
2040 ! 2. Fourier transform all fields for all latitudes
2041 !$OMP PARALLEL DO PRIVATE(DLALP,DLDALP,DLSP,DLSP2,ZFMS,ZFMA, &
2042 !$OMP INM,ILA,JM,JN,JK,JK2,JJ,JJ2,ZJM,ILONR,ILONI,FACTOR)
2043 do jm = gst(gstID)%mymBeg, gst(gstID)%mymEnd, gst(gstID)%mymSkip
2044
2045 ilonr = 2*gst(gstID)%mymIndex(jm)-1
2046 iloni = 2*gst(gstID)%mymIndex(jm)
2047 zjm = real(jm,8)
2048
2049 ! 2.1 Fetch the Legendre functions and their derivatives for this choice of "m"
2050 call getalp(dlalp,dldalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc,jm)
2051
2052 ! 2.2 Build the symmetric and anti-symmetric Fourier coefficients including
2053 ! the appropriate quadrature weights (see scientific notes)
2054 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2055 do jj = 1, gst(gstID)%njlath
2056 zfms(jj,1,jk) = 0.0d0
2057 zfms(jj,2,jk) = 0.0d0
2058 zfma(jj,1,jk) = 0.0d0
2059 zfma(jj,2,jk) = 0.0d0
2060 enddo
2061 enddo
2062
2063 ! 2.2.1 Coefficients associated with the scalar fields
2064 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2065 if(jk .gt. 2*nflev) then
2066 do jj = 1, gst(gstID)%nj
2067 jj2 = gst(gstID)%nj-jj+1
2068 if(jj.le.gst(gstID)%njlath) then
2069 ! Northern hemisphere
2070 zfms(jj,1,jk) = pgd2(ilonr,jj2,jk)
2071 zfms(jj,2,jk) = pgd2(iloni,jj2,jk)
2072 zfma(jj,1,jk) = pgd2(ilonr,jj2,jk)
2073 zfma(jj,2,jk) = pgd2(iloni,jj2,jk)
2074 else
2075 ! Southern hemisphere
2076 zfms(jj2,1,jk) = zfms(jj2,1,jk) + pgd2(ilonr,jj2,jk)
2077 zfms(jj2,2,jk) = zfms(jj2,2,jk) + pgd2(iloni,jj2,jk)
2078 zfma(jj2,1,jk) = zfma(jj2,1,jk) - pgd2(ilonr,jj2,jk)
2079 zfma(jj2,2,jk) = zfma(jj2,2,jk) - pgd2(iloni,jj2,jk)
2080 endif
2081 enddo
2082 endif
2083 enddo
2084
2085 ! 2.2.2 Coefficients associated with the wind fields
2086 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd, 2
2087 if(jk .le. 2*nflev) then
2088 jk2 = jk + 1
2089 do jj = 1, gst(gstID)%nj
2090 jj2 = gst(gstID)%nj-jj+1
2091 if(jj.le.gst(gstID)%njlath) then
2092 ! Northern hemisphere
2093 ! vorticity: symmetric coefficients
2094 zfms(jj,1,jk ) = -pgd2(iloni,jj2,jk2)
2095 zfms(jj,2,jk ) = pgd2(ilonr,jj2,jk2)
2096 ! vorticity: antisymmetric coefficients
2097 zfma(jj,1,jk ) = -pgd2(iloni,jj2,jk2)
2098 zfma(jj,2,jk ) = pgd2(ilonr,jj2,jk2)
2099 ! divergence: symmetric coefficients
2100 zfms(jj,1,jk2) = -pgd2(iloni,jj2,jk )
2101 zfms(jj,2,jk2) = pgd2(ilonr,jj2,jk )
2102 ! divergence: antisymmetric coefficients
2103 zfma(jj,1,jk2) = -pgd2(iloni,jj2,jk )
2104 zfma(jj,2,jk2) = pgd2(ilonr,jj2,jk )
2105 else
2106 ! Southern hemisphere
2107 ! vorticity: symmetric coefficients
2108 zfms(jj2,1,jk ) = zfms(jj2,1,jk ) - pgd2(iloni,jj2,jk2)
2109 zfms(jj2,2,jk ) = zfms(jj2,2,jk ) + pgd2(ilonr,jj2,jk2)
2110 ! vorticity: antisymmetric coefficients
2111 zfma(jj2,1,jk ) = zfma(jj2,1,jk ) + pgd2(iloni,jj2,jk2)
2112 zfma(jj2,2,jk ) = zfma(jj2,2,jk ) - pgd2(ilonr,jj2,jk2)
2113 ! divergence: symmetric coefficients
2114 zfms(jj2,1,jk2) = zfms(jj2,1,jk2) - pgd2(iloni,jj2,jk )
2115 zfms(jj2,2,jk2) = zfms(jj2,2,jk2) + pgd2(ilonr,jj2,jk )
2116 ! divergence: antisymmetric coefficients
2117 zfma(jj2,1,jk2) = zfma(jj2,1,jk2) + pgd2(iloni,jj2,jk )
2118 zfma(jj2,2,jk2) = zfma(jj2,2,jk2) - pgd2(ilonr,jj2,jk )
2119 endif
2120 enddo ! jj
2121 endif
2122 enddo ! jk
2123
2124 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2125 do jj = 1, gst(gstID)%njlath
2126 zfms(jj,1,jk) = dlrwt(jj)*zfms(jj,1,jk)
2127 zfms(jj,2,jk) = dlrwt(jj)*zfms(jj,2,jk)
2128 zfma(jj,1,jk) = dlrwt(jj)*zfma(jj,1,jk)
2129 zfma(jj,2,jk) = dlrwt(jj)*zfma(jj,2,jk)
2130 enddo
2131 if(jk .le. 2*nflev) then
2132 do jj = 1, gst(gstID)%njlath
2133 factor = ec_ra / cos(gst(gstID)%rlati(jj))
2134 zfms(jj,1,jk) = factor*zjm*zfms(jj,1,jk)
2135 zfms(jj,2,jk) = factor*zjm*zfms(jj,2,jk)
2136 zfma(jj,1,jk) = factor*zjm*zfma(jj,1,jk)
2137 zfma(jj,2,jk) = factor*zjm*zfma(jj,2,jk)
2138 enddo
2139 endif
2140 enddo ! jk
2141
2142 ! 2.3 First one with ALP for all scalar fields and for half the terms
2143 ! required to define the divergence and vorticity
2144 call legdir2d(jm,zfms,zfma,dlsp,dlalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc)
2145
2146 ! 2.4 Second transform with DALP to complete the construction of the
2147 ! vorticity and divergence fields
2148 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2149 if(jk .le. 2*nflev) then
2150 do jj = 1, gst(gstID)%njlath
2151 zfms(jj,1,jk) = 0.0d0
2152 zfms(jj,2,jk) = 0.0d0
2153 zfma(jj,1,jk) = 0.0d0
2154 zfma(jj,2,jk) = 0.0d0
2155 enddo
2156 endif
2157 enddo
2158
2159 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd, 2
2160 if(jk .le. 2*nflev) then
2161 jk2 = jk + 1
2162 do jj = 1, gst(gstID)%nj
2163 jj2 = gst(gstID)%nj-jj+1
2164 if(jj.le.gst(gstID)%njlath) then
2165 ! Northern hemisphere
2166 ! symmetric coefficients for zonal wind
2167 zfms(jj,1,jk ) = pgd2(ilonr,jj2,jk )
2168 zfms(jj,2,jk ) = pgd2(iloni,jj2,jk )
2169 ! antisymmetric coefficients for zonal wind
2170 zfma(jj,1,jk ) = pgd2(ilonr,jj2,jk )
2171 zfma(jj,2,jk ) = pgd2(iloni,jj2,jk )
2172 ! symmetric coefficients for meridional wind
2173 zfms(jj,1,jk2) = -pgd2(ilonr,jj2,jk2)
2174 zfms(jj,2,jk2) = -pgd2(iloni,jj2,jk2)
2175 ! antisymmetric coefficients for meridional wind
2176 zfma(jj,1,jk2) = -pgd2(ilonr,jj2,jk2)
2177 zfma(jj,2,jk2) = -pgd2(iloni,jj2,jk2)
2178 else
2179 ! Southern hemisphere
2180 ! symmetric coefficients for zonal wind
2181 zfms(jj2,1,jk ) = zfms(jj2,1,jk ) + pgd2(ilonr,jj2,jk )
2182 zfms(jj2,2,jk ) = zfms(jj2,2,jk ) + pgd2(iloni,jj2,jk )
2183 ! antisymmetric coefficients for zonal wind
2184 zfma(jj2,1,jk ) = zfma(jj2,1,jk ) - pgd2(ilonr,jj2,jk )
2185 zfma(jj2,2,jk ) = zfma(jj2,2,jk ) - pgd2(iloni,jj2,jk )
2186 ! symmetric coefficients for meridional wind
2187 zfms(jj2,1,jk2) = zfms(jj2,1,jk2) - pgd2(ilonr,jj2,jk2)
2188 zfms(jj2,2,jk2) = zfms(jj2,2,jk2) - pgd2(iloni,jj2,jk2)
2189 ! antisymmetric coefficients for meridional wind
2190 zfma(jj2,1,jk2) = zfma(jj2,1,jk2) + pgd2(ilonr,jj2,jk2)
2191 zfma(jj2,2,jk2) = zfma(jj2,2,jk2) + pgd2(iloni,jj2,jk2)
2192 endif
2193 enddo
2194 endif
2195 enddo
2196
2197 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2198 if(jk .le. 2*nflev) then
2199 do jj = 1, gst(gstID)%njlath
2200 factor = ec_ra / cos(gst(gstID)%rlati(jj))
2201 zfms(jj,1,jk) = factor*dlrwt(jj)*zfms(jj,1,jk)
2202 zfms(jj,2,jk) = factor*dlrwt(jj)*zfms(jj,2,jk)
2203 zfma(jj,1,jk) = factor*dlrwt(jj)*zfma(jj,1,jk)
2204 zfma(jj,2,jk) = factor*dlrwt(jj)*zfma(jj,2,jk)
2205 enddo
2206 endif
2207 enddo
2208
2209 call legdir2d(jm,zfma,zfms,dlsp2,dldalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc)
2210
2211 ! 2.5 Transfer the result in the global state
2212 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2213 do jn = jm, gst(gstID)%ntrunc
2214 ila = gst(gstID)%nind(jm) + jn - jm
2215 inm = jn - jm
2216 if(jk .le. 2*nflev) then
2217 psp(ila,1,jk) = -gst(gstID)%r1snp1(ila)*(dlsp(inm,1,jk) + dlsp2(inm,1,jk))
2218 psp(ila,2,jk) = -gst(gstID)%r1snp1(ila)*(dlsp(inm,2,jk) + dlsp2(inm,2,jk))
2219 else
2220 psp(ila,1,jk) = dlsp(inm,1,jk)
2221 psp(ila,2,jk) = dlsp(inm,2,jk)
2222 endif
2223 enddo
2224 enddo
2225
2226 ! End of loop on zonal wavenumbers
2227 enddo
2228 !$OMP END PARALLEL DO
2229
2230 end subroutine spgdapar
2231
2232
2233 subroutine gst_speree(psp,pgd)
2234 implicit none
2235
2236 ! Arguments:
2237 real(8), intent(in) :: psp(gst(gstID)%myNla, 2, gst(gstID)%nk)
2238 real(8), intent(out) :: pgd(gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
2239 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%nk)
2240
2241 ! Locals:
2242 real(8), allocatable :: psp2(:,:,:),pgd2(:,:,:),pgd3(:,:,:)
2243 integer :: jlat, jk, jlon
2244
2245 allocate(psp2(gst(gstID)%nla, 2, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
2246 allocate(pgd2(2*gst(gstID)%maxmcount, gst(gstID)%nj, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
2247 allocate(pgd3(gst(gstID)%ni, gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
2248
2249 ! 1.0 Transpose data along npex from N to Levels
2250 call transpose2d_NtoLev(psp,psp2)
2251
2252 ! 1.1 Inverse Legendre transform (lon -> m)
2253 call utl_tmg_start(150,'low-level--gst_lt')
2254 call spereepar(psp2,pgd2)
2255 call utl_tmg_stop(150)
2256 deallocate(psp2)
2257
2258 ! 1.2 Transpose data along npey from M to Latitudes
2259 call transpose2d_MtoLat(pgd2,pgd3)
2260 deallocate(pgd2)
2261
2262 ! 2.1 Reset to zero the modes that are not part of the truncation
2263 !$OMP PARALLEL DO PRIVATE (JLAT,JLON,JK)
2264 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2265 do jlat = gst(gstID)%myLatBeg, gst(gstID)%myLatEnd
2266 do jlon = 2*(gst(gstID)%ntrunc+1)+1, gst(gstID)%ni
2267 pgd3(jlon,jlat,jk) = 0.d0
2268 enddo
2269 enddo
2270 enddo
2271 !$OMP END PARALLEL DO
2272
2273 ! 2.2 Apply the inverse FFT
2274 call utl_tmg_start(151,'low-level--gst_fft')
2275 call fft3dvar(pgd3,+1)
2276 call utl_tmg_stop(151)
2277
2278 ! 2.3 Transpose data along npex from Levels to Longitudes
2279 call transpose2d_LevtoLon(pgd3,pgd)
2280 deallocate(pgd3)
2281
2282 end subroutine gst_speree
2283
2284
2285 subroutine gst_speree_kij(psp,pgd)
2286 implicit none
2287
2288 ! Arguments:
2289 real(8), intent(in) :: psp(gst(gstID)%myNla, 2, gst(gstID)%nk)
2290 real(8), intent(out) :: pgd(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
2291 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
2292
2293 ! Locals:
2294 real(8), allocatable :: psp2(:,:,:),pgd2(:,:,:),pgd3(:,:,:)
2295 integer :: jlat, jk, jlon, ierr
2296
2297 call utl_tmg_start(152,'low-level--gst_barr')
2298 if(mmpi_doBarrier) call rpn_comm_barrier('GRID',ierr)
2299 call utl_tmg_stop(152)
2300
2301 allocate(psp2(gst(gstID)%nla, 2, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
2302 allocate(pgd2(gst(gstID)%maxMyLevCount, 2*gst(gstID)%maxmcount, gst(gstID)%nj))
2303 allocate(pgd3(gst(gstID)%maxMyLevCount, gst(gstID)%ni, gst(gstID)%myLatBeg:gst(gstID)%myLatEnd))
2304
2305 pgd2( (gst(gstID)%myLevCount+1):gst(gstID)%maxMyLevCount, :, :) = 0.0d0
2306
2307 ! 1.0 Transpose data along npex from N to Levels
2308 call transpose2d_NtoLev(psp,psp2)
2309
2310 ! 1.1 Inverse Legendre transform (lon -> m)
2311 call utl_tmg_start(150,'low-level--gst_lt')
2312 call spereepar_kij(psp2,pgd2)
2313 call utl_tmg_stop(150)
2314 deallocate(psp2)
2315
2316 ! 1.2 Transpose data along npey from M to Latitudes
2317 call transpose2d_MtoLat_kij(pgd2,pgd3)
2318 deallocate(pgd2)
2319
2320 ! 2.1 Reset to zero the modes that are not part of the truncation
2321 !$OMP PARALLEL DO PRIVATE (JLAT,JLON,JK)
2322 do jlat = gst(gstID)%myLatBeg, gst(gstID)%myLatEnd
2323 do jlon = 2*(gst(gstID)%ntrunc+1)+1, gst(gstID)%ni
2324 do jk = 1, gst(gstID)%maxMyLevCount
2325 pgd3(jk,jlon,jlat) = 0.d0
2326 enddo
2327 enddo
2328 enddo
2329 !$OMP END PARALLEL DO
2330
2331 ! 2.2 Apply the inverse FFT
2332 call utl_tmg_start(151,'low-level--gst_fft')
2333 call fft3dvar_kij(pgd3,+1)
2334 call utl_tmg_stop(151)
2335
2336 ! 2.3 Transpose data along npex from Levels to Longitudes
2337 if( gst(gstID)%lonLatDivisible ) then
2338 if (pre_specTransReal == 4) then
2339 call transpose2d_LevtoLon_kij_mpitypes4(pgd3,pgd)
2340 else
2341 call transpose2d_LevtoLon_kij_mpitypes8(pgd3,pgd)
2342 end if
2343 else
2344 call transpose2d_LevtoLon_kij(pgd3,pgd)
2345 end if
2346 deallocate(pgd3)
2347
2348 call utl_tmg_start(152,'low-level--gst_barr')
2349 if(mmpi_doBarrier) call rpn_comm_barrier('GRID',ierr)
2350 call utl_tmg_stop(152)
2351
2352 end subroutine gst_speree_kij
2353
2354 subroutine gst_speree_ad(psp,pgd)
2355 implicit none
2356
2357 ! Arguments:
2358 real(8), intent(out) :: psp(gst(gstID)%myNla, 2, gst(gstID)%nk)
2359 real(8), intent(inout) :: pgd(gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
2360 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%nk)
2361
2362 call adjnorm(pgd)
2363
2364 call gst_reespe(psp,pgd)
2365
2366 end subroutine gst_speree_ad
2367
2368
2369 subroutine gst_reespe(psp,pgd)
2370 implicit none
2371
2372 ! Arguments:
2373 real(8), intent(out) :: psp(gst(gstID)%myNla, 2, gst(gstID)%nk)
2374 real(8), intent(in) :: pgd(gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
2375 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%nk)
2376
2377 ! Locals:
2378 real(8), allocatable :: psp2(:,:,:),pgd2(:,:,:),pgd3(:,:,:)
2379
2380 allocate(psp2(gst(gstID)%nla, 2, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
2381 allocate(pgd2(2*gst(gstID)%maxmcount, gst(gstID)%nj, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
2382 allocate(pgd3(gst(gstID)%ni, gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
2383
2384 ! Transpose data along npex from Longitudes to Levels
2385 call transpose2d_LontoLev(pgd,pgd3)
2386
2387 ! 1. Apply the FFT
2388 call utl_tmg_start(151,'low-level--gst_fft')
2389 call fft3dvar(pgd3,-1)
2390 call utl_tmg_stop(151)
2391
2392 ! Transpose data along npey from Latitudes to M
2393 call transpose2d_LattoM(pgd3,pgd2)
2394 deallocate(pgd3)
2395
2396 ! 2. Direct Legendre transform
2397 call utl_tmg_start(150,'low-level--gst_lt')
2398 call reespepar(pgd2,psp2)
2399 call utl_tmg_stop(150)
2400 deallocate(pgd2)
2401
2402 ! Transpose data along npex from Levels to N
2403 call transpose2d_LevtoN(psp2,psp)
2404 deallocate(psp2)
2405
2406 end subroutine gst_reespe
2407
2408
2409 subroutine gst_reespe_kij(psp,pgd)
2410 implicit none
2411
2412 ! Arguments:
2413 real(8), intent(out) :: psp(gst(gstID)%myNla, 2, gst(gstID)%nk)
2414 real(8), intent(in) :: pgd(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
2415 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
2416
2417 ! Locals:
2418 real(8), allocatable :: psp2(:,:,:),pgd2(:,:,:),pgd3(:,:,:)
2419 integer :: ierr
2420
2421 call utl_tmg_start(152,'low-level--gst_barr')
2422 if(mmpi_doBarrier) call rpn_comm_barrier('GRID',ierr)
2423 call utl_tmg_stop(152)
2424
2425 allocate(psp2(gst(gstID)%nla, 2, gst(gstID)%myLevBeg:gst(gstID)%myLevEnd))
2426 allocate(pgd2(gst(gstID)%maxMyLevCount, 2*gst(gstID)%maxmcount, gst(gstID)%nj))
2427 allocate(pgd3(gst(gstID)%maxMyLevCount, gst(gstID)%ni, gst(gstID)%myLatBeg:gst(gstID)%myLatEnd))
2428
2429 ! Transpose data along npex from Longitudes to Levels
2430 if( gst(gstID)%lonLatDivisible ) then
2431 if (pre_specTransReal == 4) then
2432 call transpose2d_LontoLev_kij_mpitypes4(pgd,pgd3)
2433 else
2434 call transpose2d_LontoLev_kij_mpitypes8(pgd,pgd3)
2435 end if
2436 else
2437 call transpose2d_LontoLev_kij(pgd,pgd3)
2438 end if
2439
2440 ! 1. Apply the FFT
2441 call utl_tmg_start(151,'low-level--gst_fft')
2442 call fft3dvar_kij(pgd3,-1)
2443 call utl_tmg_stop(151)
2444
2445 ! Transpose data along npey from Latitudes to M
2446 call transpose2d_LattoM_kij(pgd3,pgd2)
2447 deallocate(pgd3)
2448
2449 ! 2. Direct Legendre transform
2450 call utl_tmg_start(150,'low-level--gst_lt')
2451 call reespepar_kij(pgd2,psp2)
2452 call utl_tmg_stop(150)
2453 deallocate(pgd2)
2454
2455 ! Transpose data along npex from Levels to N
2456 call transpose2d_LevtoN(psp2,psp)
2457 deallocate(psp2)
2458
2459 call utl_tmg_start(152,'low-level--gst_barr')
2460 if(mmpi_doBarrier) call rpn_comm_barrier('GRID',ierr)
2461 call utl_tmg_stop(152)
2462
2463 end subroutine gst_reespe_kij
2464
2465 subroutine gst_speree_kij_ad(psp,pgd)
2466 implicit none
2467
2468 ! Arguments:
2469 real(8), intent(out) :: psp(gst(gstID)%myNla, 2, gst(gstID)%nk)
2470 real(8), intent(inout) :: pgd(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
2471 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
2472
2473 call adjnorm_kij(pgd)
2474
2475 call gst_reespe_kij(psp,pgd)
2476
2477 end subroutine gst_speree_kij_ad
2478
2479 subroutine adjnorm(pgd)
2480 implicit none
2481
2482 ! Arguments:
2483 real(8), intent(inout) :: pgd(gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
2484 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd, gst(gstID)%nk)
2485
2486 ! Locals:
2487 integer :: jk, jlon, jlat
2488 integer :: lat1, lat2, lon1, lon2
2489 real(8) :: rwtinv(gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
2490
2491 lat1 = gst(gstID)%myLatBeg
2492 lat2 = gst(gstID)%myLatEnd
2493 lon1 = gst(gstID)%myLonBeg
2494 lon2 = gst(gstID)%myLonEnd
2495
2496 do jlat = lat1, lat2
2497 rwtinv(jlat) = real(gst(gstID)%ni,8) / gst_getRWT(jlat, gstID)
2498 enddo
2499
2500 !$OMP PARALLEL DO PRIVATE (jk,jlat,jlon)
2501 do jk = 1, gst(gstID)%nk
2502 do jlat = lat1, lat2
2503 do jlon = lon1, lon2
2504 pgd(jlon,jlat,jk) = rwtinv(jlat) * pgd(jlon,jlat,jk)
2505 end do
2506 end do
2507 end do
2508 !$OMP END PARALLEL DO
2509
2510 end subroutine adjnorm
2511
2512 subroutine adjnorm_kij(pgd)
2513 implicit none
2514
2515 ! Arguments:
2516 real(8), intent(inout) :: pgd(gst(gstID)%nk, gst(gstID)%myLonBeg:gst(gstID)%myLonEnd, &
2517 gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
2518
2519 ! Locals:
2520 integer :: jk, jlon, jlat
2521 integer :: lat1, lat2, lon1, lon2
2522 real(8) :: rwtinv(gst(gstID)%myLatBeg:gst(gstID)%myLatEnd)
2523
2524 lat1 = gst(gstID)%myLatBeg
2525 lat2 = gst(gstID)%myLatEnd
2526 lon1 = gst(gstID)%myLonBeg
2527 lon2 = gst(gstID)%myLonEnd
2528
2529 do jlat = lat1, lat2
2530 rwtinv(jlat) = real(gst(gstID)%ni,8) / gst_getRWT(jlat, gstID)
2531 enddo
2532
2533 !$OMP PARALLEL DO PRIVATE (jlat,jlon,jk)
2534 do jlat = lat1, lat2
2535 do jlon = lon1, lon2
2536 do jk = 1, gst(gstID)%nk
2537 pgd(jk,jlon,jlat) = rwtinv(jlat) * pgd(jk,jlon,jlat)
2538 end do
2539 end do
2540 end do
2541 !$OMP END PARALLEL DO
2542
2543 end subroutine adjnorm_kij
2544
2545 subroutine spereepar(psp,pgd2)
2546 !
2547 !:Purpose: Inverse spectral transform(MPI PARALLEL LOOP)
2548 implicit none
2549
2550 ! Arguments:
2551 real(8), intent(in) :: psp(gst(gstID)%nla,2,gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2552 real(8), intent(out) :: pgd2(2*gst(gstID)%maxmcount,gst(gstID)%nj, &
2553 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2554
2555 ! Locals
2556 integer :: jj, jj2, jm, jn, ilonr, iloni, jk, ila, inm
2557 real(8) :: dlalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
2558 real(8) :: dldalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
2559 real(8) :: zfms(gst(gstID)%njlath+1,2, &
2560 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2561 real(8) :: zfma(gst(gstID)%njlath+1,2, &
2562 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2563 real(8) :: dlsp(0:gst(gstID)%ntrunc,2, &
2564 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2565
2566 ! Inverse Legendre transform
2567
2568 !$OMP PARALLEL DO PRIVATE(DLALP,DLDALP,DLSP,ZFMS,ZFMA, &
2569 !$OMP INM,ILA,JM,JN,JK,JJ,JJ2,ILONR,ILONI)
2570 do jm = gst(gstID)%mymBeg, gst(gstID)%mymEnd, gst(gstID)%mymSkip
2571
2572 ilonr = 2*gst(gstID)%mymIndex(jm)-1
2573 iloni = 2*gst(gstID)%mymIndex(jm)
2574
2575 ! 2.1 Copy global spectral state into local spectral state
2576 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2577 do jn = jm, gst(gstID)%ntrunc
2578 ila = gst(gstID)%nind(jm) + jn - jm
2579 inm = jn - jm
2580 dlsp(inm,1,jk) = psp(ila,1,jk)
2581 dlsp(inm,2,jk) = psp(ila,2,jk)
2582 enddo
2583 enddo
2584
2585 ! 2.2 Get Legendre polynomial (and its derivative) for all latitudes
2586 ! but for the chosen value of "m" from the global array
2587 call getalp(dlalp,dldalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc,jm)
2588
2589 ! 2.3 Perform the inverse Legendre transform for all fields
2590 call leginv2d(jm,zfms,zfma,dlsp,dlalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc)
2591
2592 ! 2.4 Passage to Fourier space
2593
2594 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2595 do jj = 1, gst(gstID)%nj
2596 jj2 = gst(gstID)%nj - jj + 1
2597 if(jj.le.gst(gstID)%njlath) then
2598 pgd2(ilonr,jj2,jk) = zfms(jj,1,jk) + zfma(jj,1,jk)
2599 pgd2(iloni,jj2,jk) = zfms(jj,2,jk) + zfma(jj,2,jk)
2600 else
2601 pgd2(ilonr,jj2,jk) = zfms(jj2,1,jk) - zfma(jj2,1,jk)
2602 pgd2(iloni,jj2,jk) = zfms(jj2,2,jk) - zfma(jj2,2,jk)
2603 endif
2604 enddo
2605 enddo
2606 enddo
2607 !$OMP END PARALLEL DO
2608
2609 end subroutine spereepar
2610
2611
2612 subroutine spereepar_kij(psp,pgd2)
2613 !
2614 !:Purpose: Inverse spectral transform(MPI PARALLEL LOOP)
2615 implicit none
2616
2617 ! Arguments:
2618 real(8), intent(in) :: psp(gst(gstID)%nla,2,gst(gstID)%myLevCount)
2619 real(8), intent(out) :: pgd2(gst(gstID)%maxMyLevCount,2*gst(gstID)%maxmcount, &
2620 gst(gstID)%nj)
2621
2622 ! Locals
2623 integer :: jj, jj2, jm, jn, ilonr, iloni, jk, ila, inm
2624 real(8) :: dlalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
2625 real(8) :: dldalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
2626 real(8) :: zfms(gst(gstID)%myLevCount,gst(gstID)%njlath+1,2)
2627 real(8) :: zfma(gst(gstID)%myLevCount,gst(gstID)%njlath+1,2)
2628 real(8) :: dlsp(gst(gstID)%myLevCount,0:gst(gstID)%ntrunc,2)
2629
2630 ! Inverse Legendre transform
2631
2632 !$OMP PARALLEL DO PRIVATE(DLALP,DLDALP,DLSP,ZFMS,ZFMA, &
2633 !$OMP INM,ILA,JM,JN,JK,JJ,JJ2,ILONR,ILONI)
2634 do jm = gst(gstID)%mymBeg, gst(gstID)%mymEnd, gst(gstID)%mymSkip
2635
2636 ilonr = 2*gst(gstID)%mymIndex(jm)-1
2637 iloni = 2*gst(gstID)%mymIndex(jm)
2638
2639 ! 2.1 Copy global spectral state into local spectral state
2640 do jn = jm, gst(gstID)%ntrunc
2641 ila = gst(gstID)%nind(jm) + jn - jm
2642 inm = jn - jm
2643 do jk = 1, gst(gstID)%myLevCount
2644 dlsp(jk,inm,1) = psp(ila,1,jk)
2645 dlsp(jk,inm,2) = psp(ila,2,jk)
2646 enddo
2647 enddo
2648
2649 ! 2.2 Get Legendre polynomial (and its derivative) for all latitudes
2650 ! but for the chosen value of "m" from the global array
2651 call getalp(dlalp,dldalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc,jm)
2652
2653 ! 2.3 Perform the inverse Legendre transform for all fields
2654 call leginv2d_kij(jm,zfms,zfma,dlsp,dlalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc)
2655
2656 ! 2.4 Passage to Fourier space
2657
2658 do jj = 1, gst(gstID)%nj
2659 jj2 = gst(gstID)%nj - jj + 1
2660 if(jj.le.gst(gstID)%njlath) then
2661 do jk = 1, gst(gstID)%myLevCount
2662 pgd2(jk,ilonr,jj2) = zfms(jk,jj,1) + zfma(jk,jj,1)
2663 pgd2(jk,iloni,jj2) = zfms(jk,jj,2) + zfma(jk,jj,2)
2664 enddo
2665 else
2666 do jk = 1, gst(gstID)%myLevCount
2667 pgd2(jk,ilonr,jj2) = zfms(jk,jj2,1) - zfma(jk,jj2,1)
2668 pgd2(jk,iloni,jj2) = zfms(jk,jj2,2) - zfma(jk,jj2,2)
2669 enddo
2670 endif
2671 enddo
2672
2673 enddo
2674 !$OMP END PARALLEL DO
2675
2676 end subroutine spereepar_kij
2677
2678
2679 subroutine reespepar(pgd2,psp)
2680 implicit none
2681
2682 ! Arguments:
2683 real(8), intent(out) :: psp(gst(gstID)%nla,2,gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2684 real(8), intent(in) :: pgd2(2*gst(gstID)%maxmcount,gst(gstID)%nj, &
2685 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2686
2687 ! Locals:
2688 integer :: jj,jj2,jk,ilonr, iloni
2689 integer :: jm, ila, inm, jn
2690 real(8) :: dlalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
2691 real(8) :: dldalp(0:gst(gstID)%ntrunc, gst(gstID)%njlath)
2692 real(8) :: dlsp(0:gst(gstID)%ntrunc,2, &
2693 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2694 real(8) :: zfms(gst(gstID)%njlath+1,2, &
2695 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2696 real(8) :: zfma(gst(gstID)%njlath+1,2, &
2697 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2698 real(8) :: dlrwt(gst(gstID)%nj)
2699
2700 ! 1. Adjustment needed when an odd number of latitudes is considered
2701 dlrwt(:) = gst(gstID)%rwt(:)
2702 if (mod(gst(gstID)%nj,2).ne.0) then
2703 dlrwt(gst(gstID)%njlath) = dlrwt(gst(gstID)%njlath)/2.d0
2704 end if
2705
2706 !$OMP PARALLEL DO PRIVATE(DLALP,DLDALP,DLSP,ZFMS,ZFMA, &
2707 !$OMP INM,ILA,JM,JN,JK,JJ,JJ2,ILONR,ILONI)
2708 do jm = gst(gstID)%mymBeg, gst(gstID)%mymEnd, gst(gstID)%mymSkip
2709
2710 ilonr = 2*gst(gstID)%mymIndex(jm)-1
2711 iloni = 2*gst(gstID)%mymIndex(jm)
2712
2713 ! 2.1 Fetch the Legendre functions and their derivatives for this choice of "m"
2714 call getalp(dlalp,dldalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc,jm)
2715
2716 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2717 ! 2.2 Build the symmetric and anti-symmetric Fourier coefficients including
2718 ! the appropriate quadrature weights (see scientific notes)
2719 do jj = 1, gst(gstID)%njlath
2720 zfms(jj,1,jk) = 0.0d0
2721 zfms(jj,2,jk) = 0.0d0
2722 zfma(jj,1,jk) = 0.0d0
2723 zfma(jj,2,jk) = 0.0d0
2724 enddo
2725
2726 do jj = 1, gst(gstID)%nj
2727 jj2 = gst(gstID)%nj-jj+1
2728 if(jj.le.gst(gstID)%njlath) then
2729 ! Northern hemisphere
2730 zfms(jj,1,jk) = pgd2(ilonr,jj2,jk)
2731 zfms(jj,2,jk) = pgd2(iloni,jj2,jk)
2732 zfma(jj,1,jk) = pgd2(ilonr,jj2,jk)
2733 zfma(jj,2,jk) = pgd2(iloni,jj2,jk)
2734 else
2735 ! Southern hemisphere
2736 zfms(jj2,1,jk) = zfms(jj2,1,jk) + pgd2(ilonr,jj2,jk)
2737 zfms(jj2,2,jk) = zfms(jj2,2,jk) + pgd2(iloni,jj2,jk)
2738 zfma(jj2,1,jk) = zfma(jj2,1,jk) - pgd2(ilonr,jj2,jk)
2739 zfma(jj2,2,jk) = zfma(jj2,2,jk) - pgd2(iloni,jj2,jk)
2740 endif
2741 enddo
2742
2743 do jj = 1,gst(gstID)%njlath
2744 zfms(jj,1,jk) = dlrwt(jj)*zfms(jj,1,jk)
2745 zfms(jj,2,jk) = dlrwt(jj)*zfms(jj,2,jk)
2746 zfma(jj,1,jk) = dlrwt(jj)*zfma(jj,1,jk)
2747 zfma(jj,2,jk) = dlrwt(jj)*zfma(jj,2,jk)
2748 enddo
2749
2750 enddo ! jk
2751
2752 ! 2.3 First one with ALP for all scalar fields and for half the terms
2753 ! required to define the divergence and vorticity
2754 call legdir2d(jm,zfms,zfma,dlsp,dlalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc)
2755
2756 ! 2.4 Transfer the result in the global state
2757 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2758 do jn = jm, gst(gstID)%ntrunc
2759 ila = gst(gstID)%nind(jm) + jn - jm
2760 inm = jn - jm
2761 psp(ila,1,jk) = dlsp(inm,1,jk)
2762 psp(ila,2,jk) = dlsp(inm,2,jk)
2763 enddo
2764 enddo ! jk
2765
2766 ! End of loop on zonal wavenumbers
2767 enddo
2768 !$OMP END PARALLEL DO
2769
2770 end subroutine reespepar
2771
2772
2773 subroutine reespepar_kij(pgd2,psp)
2774 implicit none
2775
2776 ! Arguments:
2777 real(8), intent(out) :: psp(gst(gstID)%nla,2,gst(gstID)%myLevCount)
2778 real(8), intent(in) :: pgd2(gst(gstID)%maxMyLevCount,2*gst(gstID)%maxmcount, &
2779 gst(gstID)%nj)
2780
2781 ! Locals:
2782 integer :: jj,jj2,jk,ilonr, iloni
2783 integer :: jm, ila, inm, jn
2784 real(8) :: dlalp(0:gst(gstID)%ntrunc,gst(gstID)%njlath)
2785 real(8) :: dldalp(0:gst(gstID)%ntrunc, gst(gstID)%njlath)
2786 real(8) :: dlsp(gst(gstID)%myLevCount,0:gst(gstID)%ntrunc,2)
2787 real(8) :: zfms(gst(gstID)%myLevCount,gst(gstID)%njlath+1,2)
2788 real(8) :: zfma(gst(gstID)%myLevCount,gst(gstID)%njlath+1,2)
2789 real(8) :: dlrwt(gst(gstID)%nj)
2790
2791 ! 1. Adjustment needed when an odd number of latitudes is considered
2792 dlrwt(:) = gst(gstID)%rwt(:)
2793 if (mod(gst(gstID)%nj,2).ne.0) then
2794 dlrwt(gst(gstID)%njlath) = dlrwt(gst(gstID)%njlath)/2.d0
2795 end if
2796
2797 !$OMP PARALLEL DO PRIVATE(DLALP,DLDALP,DLSP,ZFMS,ZFMA, &
2798 !$OMP INM,ILA,JM,JN,JK,JJ,JJ2,ILONR,ILONI)
2799 do jm = gst(gstID)%mymBeg, gst(gstID)%mymEnd, gst(gstID)%mymSkip
2800
2801 ilonr = 2*gst(gstID)%mymIndex(jm)-1
2802 iloni = 2*gst(gstID)%mymIndex(jm)
2803
2804 ! 2.1 Fetch the Legendre functions and their derivatives for this choice of "m"
2805 call getalp(dlalp,dldalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc,jm)
2806
2807 ! 2.2 Build the symmetric and anti-symmetric Fourier coefficients including
2808 ! the appropriate quadrature weights (see scientific notes)
2809 do jj = 1, gst(gstID)%njlath
2810 do jk = 1, gst(gstID)%myLevCount
2811 zfms(jk,jj,1) = 0.0d0
2812 zfms(jk,jj,2) = 0.0d0
2813 zfma(jk,jj,1) = 0.0d0
2814 zfma(jk,jj,2) = 0.0d0
2815 enddo
2816 enddo
2817
2818 do jj = 1, gst(gstID)%nj
2819 jj2 = gst(gstID)%nj-jj+1
2820 if(jj.le.gst(gstID)%njlath) then
2821 ! Northern hemisphere
2822 do jk = 1, gst(gstID)%myLevCount
2823 zfms(jk,jj,1) = pgd2(jk,ilonr,jj2)
2824 zfms(jk,jj,2) = pgd2(jk,iloni,jj2)
2825 zfma(jk,jj,1) = pgd2(jk,ilonr,jj2)
2826 zfma(jk,jj,2) = pgd2(jk,iloni,jj2)
2827 enddo
2828 else
2829 do jk = 1, gst(gstID)%myLevCount
2830 ! Southern hemisphere
2831 zfms(jk,jj2,1) = zfms(jk,jj2,1) + pgd2(jk,ilonr,jj2)
2832 zfms(jk,jj2,2) = zfms(jk,jj2,2) + pgd2(jk,iloni,jj2)
2833 zfma(jk,jj2,1) = zfma(jk,jj2,1) - pgd2(jk,ilonr,jj2)
2834 zfma(jk,jj2,2) = zfma(jk,jj2,2) - pgd2(jk,iloni,jj2)
2835 enddo
2836 endif
2837 enddo
2838
2839 do jj = 1,gst(gstID)%njlath
2840 do jk = 1, gst(gstID)%myLevCount
2841 zfms(jk,jj,1) = dlrwt(jj)*zfms(jk,jj,1)
2842 zfms(jk,jj,2) = dlrwt(jj)*zfms(jk,jj,2)
2843 zfma(jk,jj,1) = dlrwt(jj)*zfma(jk,jj,1)
2844 zfma(jk,jj,2) = dlrwt(jj)*zfma(jk,jj,2)
2845 enddo
2846 enddo
2847
2848
2849 ! 2.3 First one with ALP for all scalar fields and for half the terms
2850 ! required to define the divergence and vorticity
2851 call legdir2d_kij(jm,zfms,zfma,dlsp,dlalp,gst(gstID)%njlath,gst(gstID)%ntrunc,gst(gstID)%ntrunc)
2852
2853 ! 2.4 Transfer the result in the global state
2854 do jn = jm, gst(gstID)%ntrunc
2855 ila = gst(gstID)%nind(jm) + jn - jm
2856 inm = jn - jm
2857 do jk = 1, gst(gstID)%myLevCount
2858 psp(ila,1,jk) = dlsp(jk,inm,1)
2859 psp(ila,2,jk) = dlsp(jk,inm,2)
2860 enddo
2861 enddo
2862
2863 ! End of loop on zonal wavenumbers
2864 enddo
2865 !$OMP END PARALLEL DO
2866
2867 end subroutine reespepar_kij
2868
2869
2870 subroutine legdir2d(km,pfms,pfma,ddsp,ddalp,klath,ktrunc,ktruncdim)
2871 implicit none
2872
2873 ! Arguments:
2874 integer, intent(in) :: km
2875 integer, intent(in) :: ktrunc
2876 integer, intent(in) :: ktruncdim
2877 integer, intent(in) :: klath
2878 real(8), intent(in) :: pfms(gst(gstID)%njlath+1,2, &
2879 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2880 real(8), intent(in) :: pfma(gst(gstID)%njlath+1,2, &
2881 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2882 real(8), intent(in) :: ddalp(0:ktruncdim,klath)
2883 real(8), intent(out) :: ddsp(0:ktruncdim,2,gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2884
2885 ! Locals:
2886 integer :: jk, jlat, jn, inm, itrunc, inmp1
2887
2888 itrunc = ktrunc
2889 if(mod(ktrunc-km+1,2).eq.1) itrunc = ktrunc-1
2890
2891 if(km.ne.ktrunc)then
2892 ddsp(0:ktrunc,:,:) = 0.d0
2893 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2894 do jlat = 1, klath
2895 do jn = km, itrunc, 2
2896 inm = jn - km
2897 inmp1 = inm + 1
2898 ddsp(inm, 1,jk) = ddsp(inm, 1,jk) + ddalp(inm, jlat)*pfms(jlat,1,jk)
2899 ddsp(inmp1,1,jk) = ddsp(inmp1,1,jk) + ddalp(inmp1,jlat)*pfma(jlat,1,jk)
2900 ddsp(inm, 2,jk) = ddsp(inm, 2,jk) + ddalp(inm, jlat)*pfms(jlat,2,jk)
2901 ddsp(inmp1,2,jk) = ddsp(inmp1,2,jk) + ddalp(inmp1,jlat)*pfma(jlat,2,jk)
2902 enddo
2903 enddo
2904 enddo
2905 end if
2906
2907 if(mod(ktrunc-km+1,2).eq.1) then
2908 jn = ktrunc
2909 inm = jn - km
2910 ddsp(inm,:,:) = 0.d0
2911 do jk = gst(gstID)%myLevBeg, gst(gstID)%myLevEnd
2912 do jlat = 1, klath
2913 ddsp(inm,1,jk) = ddsp(inm,1,jk) + ddalp(inm,jlat)*pfms(jlat,1,jk )
2914 ddsp(inm,2,jk) = ddsp(inm,2,jk) + ddalp(inm,jlat)*pfms(jlat,2,jk )
2915 enddo
2916 enddo
2917 end if
2918
2919 end subroutine legdir2d
2920
2921
2922 subroutine legdir2d_kij(km,pfms,pfma,ddsp,ddalp,klath,ktrunc,ktruncdim)
2923 implicit none
2924
2925 ! Arguments:
2926 integer, intent(in) :: km
2927 integer, intent(in) :: ktrunc
2928 integer, intent(in) :: ktruncdim
2929 integer, intent(in) :: klath
2930 real(8), intent(in) :: pfms(gst(gstID)%myLevCount,gst(gstID)%njlath+1,2)
2931 real(8), intent(in) :: pfma(gst(gstID)%myLevCount,gst(gstID)%njlath+1,2)
2932 real(8), intent(in) :: ddalp(0:ktruncdim,klath)
2933 real(8), intent(out) :: ddsp(gst(gstID)%myLevCount,0:ktruncdim,2)
2934
2935 ! Locals:
2936 integer :: jk, jlat, jn, inm, itrunc, inmp1
2937
2938 itrunc = ktrunc
2939 if(mod(ktrunc-km+1,2).eq.1) itrunc = ktrunc-1
2940
2941 if(km.ne.ktrunc)then
2942 ddsp(:,0:ktrunc,:) = 0.d0
2943 do jlat = 1, klath
2944 do jn = km, itrunc, 2
2945 inm = jn - km
2946 inmp1 = inm + 1
2947 do jk = 1, gst(gstID)%myLevCount
2948 ddsp(jk,inm, 1) = ddsp(jk,inm, 1) + ddalp(inm, jlat)*pfms(jk,jlat,1)
2949 ddsp(jk,inmp1,1) = ddsp(jk,inmp1,1) + ddalp(inmp1,jlat)*pfma(jk,jlat,1)
2950 ddsp(jk,inm, 2) = ddsp(jk,inm, 2) + ddalp(inm, jlat)*pfms(jk,jlat,2)
2951 ddsp(jk,inmp1,2) = ddsp(jk,inmp1,2) + ddalp(inmp1,jlat)*pfma(jk,jlat,2)
2952 enddo
2953 enddo
2954 enddo
2955 end if
2956
2957 if(mod(ktrunc-km+1,2).eq.1) then
2958 jn = ktrunc
2959 inm = jn - km
2960 ddsp(:,inm,:) = 0.d0
2961 do jlat = 1, klath
2962 do jk = 1, gst(gstID)%myLevCount
2963 ddsp(jk,inm,1) = ddsp(jk,inm,1) + ddalp(inm,jlat)*pfms(jk,jlat,1 )
2964 ddsp(jk,inm,2) = ddsp(jk,inm,2) + ddalp(inm,jlat)*pfms(jk,jlat,2 )
2965 enddo
2966 enddo
2967 end if
2968
2969 end subroutine legdir2d_kij
2970
2971
2972 subroutine leginv2d(km,pfms,pfma,ddsp,ddalp,klath,ktrunc,ktruncdim)
2973 implicit none
2974
2975 ! Arguments:
2976 integer, intent(in) :: km
2977 integer, intent(in) :: ktrunc
2978 integer, intent(in) :: ktruncdim
2979 integer, intent(in) :: klath
2980 real(8), intent(out) :: pfms(gst(gstID)%njlath+1,2, &
2981 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2982 real(8), intent(out) :: pfma(gst(gstID)%njlath+1,2, &
2983 gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2984 real(8), intent(in) :: ddalp(0:ktruncdim,klath)
2985 real(8), intent(in) :: ddsp(0:ktruncdim,2,gst(gstID)%myLevBeg:gst(gstID)%myLevEnd)
2986
2987 ! Locals:
2988 integer :: jk, jlat, jn, inm, itrunc, inmp1
2989
2990 itrunc = ktrunc
2991 if(mod(ktrunc-km+1,2).eq.1) itrunc = ktrunc-1
2992
2993 if(km .ne. ktrunc)then
2994 do jk = gst(gstID)%myLevBeg,gst(gstID)%myLevEnd
2995 pfms(:,:,jk) = 0.d0
2996 pfma(:,:,jk) = 0.d0
2997 do jn = km, itrunc, 2
2998 inm = jn - km
2999 inmp1 = inm + 1
3000 do jlat = 1, klath
3001 pfms(jlat,1,jk) = pfms(jlat,1,jk) + ddalp(inm, jlat) * ddsp(inm, 1,jk)
3002 pfma(jlat,1,jk) = pfma(jlat,1,jk) + ddalp(inmp1,jlat) * ddsp(inmp1,1,jk)
3003 pfms(jlat,2,jk) = pfms(jlat,2,jk) + ddalp(inm, jlat) * ddsp(inm, 2,jk)
3004 pfma(jlat,2,jk) = pfma(jlat,2,jk) + ddalp(inmp1,jlat) * ddsp(inmp1,2,jk)
3005 enddo
3006 enddo
3007 enddo
3008 end if
3009
3010 if(mod(ktrunc-km+1,2).eq.1) then
3011 jn = ktrunc
3012 if (km .ne. ktrunc) then
3013 inm = jn - km
3014 do jk = gst(gstID)%myLevBeg,gst(gstID)%myLevEnd
3015 do jlat = 1, klath
3016 pfms(jlat,1,jk) = pfms(jlat,1,jk) + ddalp(inm,jlat) * ddsp(inm,1,jk)
3017 pfms(jlat,2,jk) = pfms(jlat,2,jk) + ddalp(inm,jlat) * ddsp(inm,2,jk)
3018 enddo
3019 enddo
3020 else
3021 inm = jn - km
3022 do jk = gst(gstID)%myLevBeg,gst(gstID)%myLevEnd
3023 pfms(:,:,jk) = 0.d0
3024 pfma(:,:,jk) = 0.d0
3025 do jlat = 1, klath
3026 pfms(jlat,1,jk) = ddalp(inm,jlat) * ddsp(inm,1,jk)
3027 pfms(jlat,2,jk) = ddalp(inm,jlat) * ddsp(inm,2,jk)
3028 enddo
3029 enddo
3030 end if
3031 end if
3032
3033 end subroutine leginv2d
3034
3035
3036 subroutine leginv2d_kij(km,pfms,pfma,ddsp,ddalp,klath,ktrunc,ktruncdim)
3037 implicit none
3038
3039 ! Arguments:
3040 integer, intent(in) :: km
3041 integer, intent(in) :: ktrunc
3042 integer, intent(in) :: ktruncdim
3043 integer, intent(in) :: klath
3044 real(8), intent(out) :: pfms(gst(gstID)%myLevCount,gst(gstID)%njlath+1,2)
3045 real(8), intent(out) :: pfma(gst(gstID)%myLevCount,gst(gstID)%njlath+1,2)
3046 real(8), intent(in) :: ddalp(0:ktruncdim,klath)
3047 real(8), intent(in) :: ddsp(gst(gstID)%myLevCount,0:ktruncdim,2)
3048
3049 ! Locals:
3050 integer :: jk, jlat, jn, inm, itrunc, inmp1
3051
3052 itrunc = ktrunc
3053 if(mod(ktrunc-km+1,2).eq.1) itrunc = ktrunc-1
3054
3055 if(km .ne. ktrunc)then
3056 pfms(:,:,:) = 0.d0
3057 pfma(:,:,:) = 0.d0
3058 do jn = km, itrunc, 2
3059 inm = jn - km
3060 inmp1 = inm + 1
3061 do jlat = 1, klath
3062 do jk = 1, gst(gstID)%myLevCount
3063 pfms(jk,jlat,1) = pfms(jk,jlat,1) + ddalp(inm, jlat) * ddsp(jk,inm, 1)
3064 pfma(jk,jlat,1) = pfma(jk,jlat,1) + ddalp(inmp1,jlat) * ddsp(jk,inmp1,1)
3065 pfms(jk,jlat,2) = pfms(jk,jlat,2) + ddalp(inm, jlat) * ddsp(jk,inm, 2)
3066 pfma(jk,jlat,2) = pfma(jk,jlat,2) + ddalp(inmp1,jlat) * ddsp(jk,inmp1,2)
3067 enddo
3068 enddo
3069 enddo
3070 end if
3071
3072 if(mod(ktrunc-km+1,2).eq.1) then
3073 jn = ktrunc
3074 if (km .ne. ktrunc) then
3075 inm = jn - km
3076 do jlat = 1, klath
3077 do jk = 1, gst(gstID)%myLevCount
3078 pfms(jk,jlat,1) = pfms(jk,jlat,1) + ddalp(inm,jlat) * ddsp(jk,inm,1)
3079 pfms(jk,jlat,2) = pfms(jk,jlat,2) + ddalp(inm,jlat) * ddsp(jk,inm,2)
3080 enddo
3081 enddo
3082 else
3083 inm = jn - km
3084 pfms(:,:,:) = 0.d0
3085 pfma(:,:,:) = 0.d0
3086 do jlat = 1, klath
3087 do jk = 1, gst(gstID)%myLevCount
3088 pfms(jk,jlat,1) = ddalp(inm,jlat) * ddsp(jk,inm,1)
3089 pfms(jk,jlat,2) = ddalp(inm,jlat) * ddsp(jk,inm,2)
3090 enddo
3091 enddo
3092 end if
3093 end if
3094
3095 end subroutine leginv2d_kij
3096
3097
3098 subroutine allocate_comleg
3099 !
3100 !:Purpose: Subroutine for initializing the Legendre transform
3101 implicit none
3102
3103 allocate(gst(gstID)%rmu(gst(gstID)%nj))
3104 allocate(gst(gstID)%rwt(gst(gstID)%nj))
3105 allocate(gst(gstID)%rwocs(gst(gstID)%nj))
3106 allocate(gst(gstID)%r1mu2(gst(gstID)%nj))
3107 allocate(gst(gstID)%rsqm2(gst(gstID)%nj))
3108 allocate(gst(gstID)%rcolat(gst(gstID)%nj))
3109 allocate(gst(gstID)%r1qm2(gst(gstID)%nj))
3110 allocate(gst(gstID)%r1mui(gst(gstID)%nj))
3111 allocate(gst(gstID)%r1mua(gst(gstID)%nj))
3112 allocate(gst(gstID)%rlati((-1):(gst(gstID)%nj+2)))
3113 allocate(gst(gstID)%nind(0:gst(gstID)%ntrunc))
3114 allocate(gst(gstID)%nindrh(0:gst(gstID)%ntrunc))
3115 allocate(gst(gstID)%nclm(0:gst(gstID)%ntrunc))
3116
3117 end subroutine allocate_comleg
3118
3119
3120 subroutine suleg(lverbose_opt)
3121 !
3122 !:Purpose: To initializethe Gaussian latitudes, weights and related
3123 ! quantities
3124 implicit none
3125
3126 ! Arguments:
3127 logical, optional, intent(in) :: lverbose_opt
3128
3129 ! Locals:
3130 logical :: lverbose
3131 integer :: jlat, jm
3132 real(8) :: zpisu2
3133 external gauss
3134
3135 ! Some explanation:
3136 ! rmu = sin(latitude)
3137 ! colat = pi/2 - latitude
3138 ! rwt = gauss weight (complicated)
3139 ! rwocs = rwt / cos(latitude)^2
3140
3141 if(present(lverbose_opt)) then
3142 lverbose = lverbose_opt
3143 else
3144 lverbose = .false.
3145 endif
3146
3147 if(mmpi_myid.eq.0) write(*,fmt='(//,6(" ***********"))')
3148 if(mmpi_myid.eq.0) write(*,*)' SULEG: initialisation of Gaussian', &
3149 ' latitudes, weights, etc...'
3150 if(mmpi_myid.eq.0) write(*,fmt='(6(" ***********"))')
3151
3152 ! 1. GAUSSIAN LATITUDES AND WEIGHTS OVER AN HEMISPHERE
3153 ! -------------------------------------------------
3154
3155 call gauss8(gst(gstID)%njlath,gst(gstID)%rmu(1),gst(gstID)%rwt(1), &
3156 gst(gstID)%rsqm2(1),gst(gstID)%rcolat(1),gst(gstID)%rwocs(1), &
3157 gst(gstID)%r1qm2(1),gst(gstID)%r1mui(1),gst(gstID)%r1mu2(1))
3158
3159 do jlat = 1, gst(gstID)%njlath
3160 gst(gstID)%rlati(jlat) = asin(gst(gstID)%rmu(jlat))
3161 gst(gstID)%r1mua(jlat) = ec_r1sa*gst(gstID)%r1mui(jlat)
3162 enddo
3163
3164 ! 2. COMPLETION FOR THE SOUTHERN HEMISPHERE
3165 ! --------------------------------------
3166
3167 do jlat = gst(gstID)%njlath +1, gst(gstID)%nj
3168 gst(gstID)%rmu(jlat) = -gst(gstID)%rmu(2*gst(gstID)%njlath +1 - jlat)
3169 gst(gstID)%rwocs(jlat) = gst(gstID)%rwocs(2*gst(gstID)%njlath +1 - jlat)
3170 gst(gstID)%r1mu2(jlat) = gst(gstID)%r1mu2(2*gst(gstID)%njlath +1 - jlat)
3171 gst(gstID)%rsqm2(jlat) = gst(gstID)%rsqm2(2*gst(gstID)%njlath +1 - jlat)
3172 gst(gstID)%r1qm2(jlat) = gst(gstID)%r1qm2(2*gst(gstID)%njlath +1 - jlat)
3173 gst(gstID)%r1mui(jlat) = gst(gstID)%r1mui(2*gst(gstID)%njlath +1 - jlat)
3174 gst(gstID)%r1mua(jlat) = gst(gstID)%r1mua(2*gst(gstID)%njlath +1 - jlat)
3175 gst(gstID)%rwt(jlat) = gst(gstID)%rwt(2*gst(gstID)%njlath +1 - jlat)
3176 gst(gstID)%rlati(jlat) = - gst(gstID)%rlati (2*gst(gstID)%njlath +1 - jlat)
3177 enddo
3178
3179 zpisu2 = MPC_PI_R8/2.d0
3180 do jlat = 1, gst(gstID)%nj
3181 gst(gstID)%rcolat(jlat) = zpisu2 - gst(gstID)%rlati(jlat)
3182 enddo
3183
3184 !* 3. Overdimensioning for interpolation
3185
3186 gst(gstID)%rlati(-1) = MPC_PI_R8-gst(gstID)%rlati(1)
3187 gst(gstID)%rlati(0) = MPC_PI_R8*.5d0
3188 gst(gstID)%rlati(gst(gstID)%nj+1) =-MPC_PI_R8*.5d0
3189 gst(gstID)%rlati(gst(gstID)%nj+2) =-MPC_PI_R8-gst(gstID)%rlati(gst(gstID)%nj)
3190
3191 !* 4. Print the content of GAUS
3192
3193 if(lverbose.and.mmpi_myid.eq.0) write(*,fmt='(" JLAT:",4X," RLATI",8X,"RCOLAT",8X,"RMU",10X ,"RWT",12X,"RW0CS")')
3194 do jlat = 1, gst(gstID)%nj
3195 if(lverbose.and.mmpi_myid.eq.0) write(*,fmt='(2X,I4,5(2X,G23.16))') &
3196 jlat,gst(gstID)%rlati(jlat),gst(gstID)%rcolat(jlat), gst(gstID)%rmu(jlat) &
3197 ,gst(gstID)%rwt(jlat),gst(gstID)%rwocs(jlat)
3198 enddo
3199
3200 if(lverbose.and.mmpi_myid.eq.0) write(*,fmt='(//," JLAT:",4X,"R1MU2",8X,"RSQM2",9X,"R1QM2",10X,"R1MUI",10X,"R1MUA")')
3201
3202 do jlat = 1, gst(gstID)%nj
3203 if(lverbose.and.mmpi_myid.eq.0) &
3204 write(*,fmt='(2X,I4,5(2X,G23.16))') jlat,gst(gstID)%r1mu2(jlat),gst(gstID)%rsqm2(jlat),gst(gstID)%r1qm2(jlat) &
3205 ,gst(gstID)%r1mui(jlat),gst(gstID)%r1mua(jlat)
3206 enddo
3207
3208 !* 5. Positioning within spectral arrays
3209
3210 do jm = 0, gst(gstID)%ntrunc
3211 gst(gstID)%nind(jm) = jm*(gst(gstID)%ntrunc+1) - (jm*(jm-1))/2 + 1
3212 gst(gstID)%nindrh(jm) = jm*(gst(gstID)%ntrunc+1) + 1
3213 gst(gstID)%nclm(jm) = gst(gstID)%ntrunc - jm + 1
3214 enddo
3215
3216 if(lverbose.and.mmpi_myid.eq.0) write(*,fmt='(/," NIND(0:NTRUNC):",/,10(2X,I8))') &
3217 (gst(gstID)%nind(jm),jm=0,gst(gstID)%ntrunc)
3218 if(lverbose.and.mmpi_myid.eq.0) write(*,fmt='(" NINDRH(0:NTRUNC):",/,10(2X,I8))') &
3219 (gst(gstID)%nindrh(jm),jm=0,gst(gstID)%ntrunc)
3220 if(lverbose.and.mmpi_myid.eq.0) write(*,fmt='(" NCLM(0:NTRUNC):",/,10(2X,I8))') &
3221 (gst(gstID)%nclm(jm),jm=0,gst(gstID)%ntrunc)
3222
3223 end subroutine suleg
3224
3225
3226 subroutine gauss8(nracp,racp,pg,sia,rad,pgssin2,sinm1,sinm2,sin2)
3227 implicit none
3228
3229 ! Arguments:
3230 integer, intent(in) :: nracp
3231 real(8), intent(out) :: racp(*)
3232 real(8), intent(out) :: pg(*)
3233 real(8), intent(out) :: sia(*)
3234 real(8), intent(out) :: rad(*)
3235 real(8), intent(out) :: pgssin2(*)
3236 real(8), intent(out) :: sinm1(*)
3237 real(8), intent(out) :: sinm2(*)
3238 real(8), intent(out) :: sin2(*)
3239
3240 ! Locals:
3241 real(8) :: xlim,pi,fi,fi1,fn,dot,dn,dn1,a,b,c,g,gm,gp,gt,ractemp,gtemp
3242 integer :: i,ir,irm,irp
3243
3244 xlim = 1.d-13
3245 pi = 4.d0*atan(1.d0)
3246 ir = 2*nracp
3247 fi = dble(ir)
3248 fi1 = fi+1.d0
3249 fn = dble(nracp)
3250
3251 do i = 1,nracp
3252 dot = dble(i-1)
3253 racp(i) = -pi*.5d0*(dot+.5d0)/fn + pi*.5d0
3254 racp(i) = sin(racp(i))
3255 enddo
3256
3257 dn = fi/sqrt(4.d0*fi*fi-1.d0)
3258 dn1 = fi1/sqrt(4.d0*fi1*fi1-1.d0)
3259 a = dn1*fi
3260 b = dn*fi1
3261 irp = ir + 1
3262 irm = ir -1
3263
3264 do i = 1,nracp
326542 call ordleg8(g,racp(i),ir)
3266 call ordleg8(gm,racp(i),irm)
3267 call ordleg8(gp,racp(i),irp)
3268 gt = (a*gp-b*gm)/(racp(i)*racp(i)-1.d0)
3269 ractemp = racp(i) - g/gt
3270 gtemp = racp(i) - ractemp
3271 racp(i) = ractemp
3272 if( abs(gtemp).gt.xlim) go to 42
3273 enddo
3274
3275 do i = 1,nracp
3276 a = 2.d0*(1.d0-racp(i)**2)
3277 call ordleg8(b,racp(i),irm)
3278 b = b*b*fi*fi
3279 pg(i) = a*(fi-.5d0)/b
3280 rad(i) = acos(racp(i))
3281 sia(i) = sin(rad(i))
3282 c = (sia(i))**2
3283 sinm1(i) = 1.d0/sia(i)
3284 sinm2(i) = 1.d0/c
3285 pgssin2(i) = pg(i)/c
3286 sin2(i) = c
3287 enddo
3288
3289 end subroutine gauss8
3290
3291
3292 subroutine ordleg8(sx,coa,ir)
3293 implicit none
3294
3295 ! Arguments:
3296 real(8), intent(out) :: sx
3297 real(8), intent(in) :: coa
3298 integer, intent(in) :: ir
3299
3300 ! Locals:
3301 integer :: n,kk,k,n1,irpp,irppm
3302 real(8) :: pi,sqr2,delta,sia,theta,c1,c4,s1,ang,fk,fn,fn2,fn2sq,a,b
3303
3304 pi = 4.d0*atan(1.d0)
3305 sqr2 = sqrt(2.d0)
3306 irpp = ir + 1
3307 irppm = irpp - 1
3308 delta = acos(coa)
3309 sia = sin(delta)
3310
3311 theta = delta
3312 c1 = sqr2
3313
3314 do n = 1,irppm
3315 fn2 = dble(2*n)
3316 fn2sq = fn2*fn2
3317 c1 = c1*sqrt(1.d0 - 1.d0/fn2sq)
3318 enddo
3319
3320 n = irppm
3321 fn = dble(n)
3322 ang = fn*theta
3323 s1 = 0.d0
3324 c4 = 1.d0
3325 a =-1.d0
3326 b = 0.d0
3327 n1 = n+1
3328
3329 do kk = 1,n1,2
3330 k = kk-1
3331 if (k.eq.n) c4 = 0.5d0*c4
3332 s1 = s1+c4* cos(ang)
3333 a = a+2.d0
3334 b = b+1.d0
3335 fk = dble(k)
3336 ang = theta*(fn-fk-2.d0)
3337 c4 = ( a * (fn-b+1.d0) / (b*(fn2-a)) )*c4
3338 enddo
3339
3340 sx = s1*c1
3341
3342 end subroutine ordleg8
3343
3344
3345 subroutine sualp
3346 implicit none
3347
3348 ! Locals:
3349 integer :: jj,jlat,jm,jn,ilat
3350 integer :: ilarh,ila,ilatbd
3351 real(8) :: dlalp(gst(gstID)%nlarh,nlatbd)
3352 real(8) :: dldalp(gst(gstID)%nlarh,nlatbd)
3353 !
3354 ! Memory allocation for Legendre polynomials
3355 !
3356 if(mmpi_myid.eq.0) write(*,*) 'allocating dalp:',gst(gstID)%nla,gst(gstID)%njlath,gst(gstID)%nla*gst(gstID)%njlath
3357 allocate( gst(gstID)%dalp(gst(gstID)%nla,gst(gstID)%njlath) )
3358 allocate( gst(gstID)%dealp(gst(gstID)%nla,gst(gstID)%njlath))
3359 if(mmpi_myid.eq.0) write(*,*) 'succeeded'
3360
3361 latitudes: do jlat = 1, gst(gstID)%njlath, nlatbd
3362 ilatbd = min(nlatbd,gst(gstID)%njlath - jlat + 1)
3363
3364 if(ilatbd.eq.8) then
3365 call allp(dlalp,dldalp,gst(gstID)%rmu(jlat),gst(gstID)%nclm(0),gst(gstID)%ntrunc,ilatbd)
3366 else
3367 call allp2(dlalp,dldalp,gst(gstID)%rmu(jlat),gst(gstID)%ntrunc,ilatbd)
3368 endif
3369
3370 do jm = 0,gst(gstID)%ntrunc
3371 do jn = jm,gst(gstID)%ntrunc
3372 ila = gst(gstID)%nind(jm) + jn -jm
3373 ilarh = gst(gstID)%nindrh(jm) + jn-jm
3374 do jj = 1,ilatbd
3375 ilat = jlat+jj-1
3376 gst(gstID)%dalp (ila,jlat+jj-1) = dlalp (ilarh,jj)
3377 gst(gstID)%dealp(ila,jlat+jj-1) = dldalp(ilarh,jj)
3378 enddo
3379 enddo
3380 enddo
3381 enddo latitudes
3382
3383 end subroutine sualp
3384
3385
3386 subroutine getalp(ddalp,dddalp,klath,ktrunc,ktruncdim ,km)
3387 implicit none
3388
3389 ! Arguments:
3390 real(8) :: ddalp(0:ktruncdim,klath)
3391 real(8) :: dddalp(0:ktruncdim,klath)
3392 integer :: klath
3393 integer :: ktrunc
3394 integer :: ktruncdim
3395 integer :: km
3396
3397 ! Locals:
3398 integer :: ila,ind
3399 integer :: jlat,jn, jlen
3400
3401 do jlat = 1,klath
3402 do jlen = 0, ktrunc
3403 ddalp(jlen,jlat) = 0.d0
3404 dddalp(jlen,jlat)= 0.d0
3405 end do
3406 end do
3407
3408 do jlat = 1, klath
3409 do jn = km, ktrunc
3410 ila = gst(gstID)%nind(km) + jn-km
3411 ind = jn-km
3412 ddalp(ind,jlat) = gst(gstID)%dalp(ila,jlat)
3413 dddalp(ind,jlat) = gst(gstID)%dealp(ila,jlat)
3414 end do
3415 end do
3416
3417 end subroutine getalp
3418
3419
3420 subroutine allp( p , g , x , lr , r , nlatp)
3421
3422 implicit none
3423
3424 ! Arguments:
3425 integer :: r
3426 integer :: nlatp
3427 integer :: lr(0:r)
3428 real(8) :: p(0:r,0:r,nlatp)
3429 real(8) :: g(0:r,0:r,nlatp)
3430 real(8) :: x(nlatp)
3431
3432 ! Locals:
3433 real(8) :: onehalf
3434 real(8) :: xp , xp2, p0, enm, fnm
3435 integer :: ilat , m , l , n
3436
3437 data onehalf /0.5d0/
3438
3439 do ilat = 1,nlatp
3440 xp2 = sqrt( 1.0d0 - x(ilat) ** 2 )
3441 p(0,0,ilat) = sqrt(onehalf)
3442 do m = 1,r
3443 xp = real(m,8)
3444 p(0,m,ilat) = sqrt( (2.0d0*xp+1.0d0)/(2.0d0*xp) ) * xp2 * p(0,m-1,ilat)
3445 enddo
3446 enddo
3447
3448 do ilat = 1,nlatp
3449 do m = 0,r
3450 xp = real(m,8)
3451 g(0,m,ilat) = - x(ilat)*xp * p(0,m,ilat)
3452 enddo
3453 enddo
3454 do n = 1,r
3455 do m = 0,lr(n)-1
3456 l = 1
3457 p0 = real(m+n,8)
3458 xp = real(m,8)
3459 enm = sqrt( ((p0*p0-xp*xp)*(2.0d0*p0+1.0d0))/(2.0d0*p0-1.0d0) )
3460 fnm = sqrt( (2.0d0*p0+1.0d0)/((p0*p0-xp*xp)*(2.0d0*p0-1.0d0)) )
3461
3462 p(n,m,l) = ( x(l) * p0 * p(n-1,m,l) - g(n-1,m,l) ) * fnm
3463 g(n,m,l) = enm * p(n-1,m,l) - x(l) * p0 * p(n,m,l)
3464 l = l + 1
3465 p(n,m,l) = ( x(l) * p0 * p(n-1,m,l) - g(n-1,m,l) ) * fnm
3466 g(n,m,l) = enm * p(n-1,m,l) - x(l) * p0 * p(n,m,l)
3467 l = l + 1
3468 p(n,m,l) = ( x(l) * p0 * p(n-1,m,l) - g(n-1,m,l) ) * fnm
3469 g(n,m,l) = enm * p(n-1,m,l) - x(l) * p0 * p(n,m,l)
3470 l = l + 1
3471 p(n,m,l) = ( x(l) * p0 * p(n-1,m,l) - g(n-1,m,l) ) * fnm
3472 g(n,m,l) = enm * p(n-1,m,l) - x(l) * p0 * p(n,m,l)
3473 l = l + 1
3474 p(n,m,l) = ( x(l) * p0 * p(n-1,m,l) - g(n-1,m,l) ) * fnm
3475 g(n,m,l) = enm * p(n-1,m,l) - x(l) * p0 * p(n,m,l)
3476 l = l + 1
3477 p(n,m,l) = ( x(l) * p0 * p(n-1,m,l) - g(n-1,m,l) ) * fnm
3478 g(n,m,l) = enm * p(n-1,m,l) - x(l) * p0 * p(n,m,l)
3479 l = l + 1
3480 p(n,m,l) = ( x(l) * p0 * p(n-1,m,l) - g(n-1,m,l) ) * fnm
3481 g(n,m,l) = enm * p(n-1,m,l) - x(l) * p0 * p(n,m,l)
3482 l = l + 1
3483 p(n,m,l) = ( x(l) * p0 * p(n-1,m,l) - g(n-1,m,l) ) * fnm
3484 g(n,m,l) = enm * p(n-1,m,l) - x(l) * p0 * p(n,m,l)
3485 l = l + 1
3486 enddo
3487 enddo
3488
3489 end subroutine allp
3490
3491
3492 subroutine allp2( p , g , x , r , nlatp)
3493
3494 implicit none
3495
3496 ! Arguments:
3497 integer :: r
3498 integer :: nlatp
3499 real(8) :: p(0:r,0:r,nlatp)
3500 real(8) :: g(0:r,0:r,nlatp)
3501 real(8) :: x(nlatp)
3502
3503 ! Locals:
3504 real(8) :: onehalf
3505 real(8) :: xp , xp2, p0, enm, fnm
3506 integer :: ilat , m , l , n, jlat
3507
3508 data onehalf /0.5d0/
3509
3510 do ilat = 1,nlatp
3511 xp2 = sqrt( 1.0d0 - x(ilat) ** 2 )
3512 p(0,0,ilat) = sqrt(onehalf)
3513 do m = 1,r
3514 xp = real(m,8)
3515 p(0,m,ilat) = sqrt( (2.0d0*xp+1.0d0)/(2.0d0*xp) ) * xp2 * p(0,m-1,ilat)
3516 enddo
3517 enddo
3518
3519 do ilat = 1,nlatp
3520 do m = 0,r
3521 xp = real(m,8)
3522 g(0,m,ilat) = - x(ilat)*xp * p(0,m,ilat)
3523 enddo
3524 enddo
3525
3526 do n = 1,r
3527 do m = 0, r
3528 p0 = real(m+n,8)
3529 xp = real(m,8)
3530 enm = sqrt( ((p0*p0-xp*xp)*(2.0d0*p0+1.0d0))/(2.0d0*p0-1.0d0) )
3531 fnm = sqrt( (2.0d0*p0+1.0d0)/((p0*p0-xp*xp)*(2.0d0*p0-1.0d0)) )
3532
3533 do jlat = 1, nlatp
3534 l = jlat
3535 p(n,m,l) = ( x(l) * p0 * p(n-1,m,l) - g(n-1,m,l) ) * fnm
3536 g(n,m,l) = enm * p(n-1,m,l) - x(l) * p0 * p(n,m,l)
3537 enddo
3538 enddo
3539 enddo
3540
3541 end subroutine allp2
3542
3543
3544 subroutine gst_zlegpol(gstID_in)
3545 !
3546 !:Purpose: To evaluate Legendre polynomials restricted to (n,m) = (n,0)
3547 implicit none
3548
3549 ! Arguments:
3550 integer :: gstID_in
3551
3552 ! Loclas:
3553 integer :: jn, jlat
3554 real(8) :: dlfact1, dlfact2, dln
3555 real(8) :: dlnorm(0:gst(gstID)%ntrunc)
3556
3557 allocate(gst(gstID_in)%zleg(0:gst(gstID)%ntrunc,gst(gstID)%nj))
3558
3559 do jlat = 1, gst(gstID_in)%nj
3560 gst(gstID_in)%zleg(0,gst(gstID_in)%nj-jlat+1) = sqrt(0.5d0)
3561 gst(gstID_in)%zleg(1,gst(gstID_in)%nj-jlat+1) = sqrt(1.5d0)*gst(gstID_in)%rmu(jlat)
3562 enddo
3563
3564 do jn = 0, gst(gstID_in)%ntrunc
3565 dln = 1.d0*real(jn,8)
3566 dlnorm(jn) = dsqrt((2.d0*dln + 1.d0)/2.d0)
3567 enddo
3568
3569 do jn = 1, gst(gstID_in)%ntrunc-1
3570 dln = real(jn,8)
3571 dlfact1 = ((2.d0*dln+1.d0)/(dln+1.d0))*(dlnorm(jn+1)/dlnorm(jn))
3572 dlfact2 = (dln/(dln+1.d0))*(dlnorm(jn+1)/dlnorm(jn-1))
3573 do jlat = 1,gst(gstID_in)%nj
3574 gst(gstID_in)%zleg(jn+1,gst(gstID_in)%nj-jlat+1) = &
3575 dlfact1*gst(gstID_in)%rmu(jlat)*dble(gst(gstID_in)%zleg(jn,gst(gstID_in)%nj-jlat+1)) &
3576 - dlfact2*dble(gst(gstID_in)%zleg(jn-1,gst(gstID_in)%nj-jlat+1))
3577 enddo
3578 enddo
3579
3580 end subroutine gst_zlegpol
3581
3582
3583 subroutine gst_zlegdir(gstID_in,pf,pn,klev)
3584 !
3585 !:Purpose: Direct Legendre transform restricted to
3586 !
3587 implicit none
3588
3589 ! Arguments:
3590 integer :: gstID_in
3591 real(8) :: pf(gst(gstID_in)%nj,klev) ! PF(NJ,KLEV): field in physical space
3592 real(8) :: pn(0:gst(gstID_in)%ntrunc,klev) ! PN(0:ntrunc, KLEV): spectral coefficients
3593 integer :: klev ! number of fields to transform
3594
3595 ! Locals:
3596 integer :: jlat, jn
3597 real(8), allocatable :: zwork(:,:)
3598
3599 allocate(zwork(0:gst(gstID_in)%ntrunc,gst(gstID_in)%nj))
3600
3601 do jlat = 1, gst(gstID_in)%nj
3602 do jn = 0, gst(gstID_in)%ntrunc
3603 zwork(jn,jlat) = gst(gstID_in)%zleg(jn,jlat)*gst_getRWT(jlat,gstID_in)
3604 end do
3605 end do
3606
3607 call dgemm('N','N',gst(gstID_in)%ntrunc+1, klev, gst(gstID_in)%nj, 1.0d0, zwork(0,1), &
3608 gst(gstID_in)%ntrunc+1, pf(1,1), &
3609 gst(gstID_in)%nj, 0.0d0, pn(0,1), gst(gstID_in)%ntrunc+1)
3610
3611 deallocate(zwork)
3612
3613 end subroutine gst_zlegdir
3614
3615
3616 subroutine gst_zleginv(gstID_in,pf,pn,klev)
3617 !
3618 !:Purpose: Direct Legendre transform restricted to fields that vary with
3619 ! latitude only
3620 !
3621 implicit none
3622
3623 ! Arguments:
3624 integer :: gstID_in
3625 real(8) :: pf(gst(gstID_in)%nj,klev) ! PF(KNJDIM,KLEVDIM) : field in physical space
3626 real(8) :: pn(0:gst(gstID_in)%ntrunc,klev) ! PN(0:KNDIM, KLEVDIM): spectral coefficients
3627 integer :: klev ! number of fields to transform
3628
3629 ! Locals:
3630 integer :: jlat, jn
3631 real(8), allocatable :: zwork(:,:)
3632
3633 allocate(zwork(0:gst(gstID_in)%ntrunc,gst(gstID_in)%nj))
3634
3635 do jlat = 1, gst(gstID_in)%nj
3636 do jn = 0, gst(gstID_in)%ntrunc
3637 zwork(jn,jlat) = gst(gstID_in)%zleg(jn,jlat)
3638 end do
3639 end do
3640
3641 call dgemm('T','N',gst(gstID_in)%nj, klev, gst(gstID_in)%ntrunc+1, 1.0d0, zwork(0,1), &
3642 gst(gstID_in)%ntrunc+1, pn(0,1), &
3643 gst(gstID_in)%ntrunc+1, 0.0d0, pf(1,1), gst(gstID_in)%nj)
3644
3645 deallocate(zwork)
3646
3647 end subroutine gst_zleginv
3648
3649
3650! ---------------------------------------------
3651! FFT subroutines
3652! ---------------------------------------------
3653
3654 subroutine fft3dvar(pgd,kdir)
3655 implicit none
3656
3657 ! Arguments:
3658 real(8) :: pgd(:,:,:)
3659 integer :: kdir
3660
3661 ! Locals:
3662 integer :: kfield,ni,nj
3663 real(8), allocatable :: pgd2(:,:)
3664 integer :: ji,jj,jk,ijump,i
3665
3666 ni = size(pgd,1)
3667 nj = size(pgd,2)
3668 kfield = size(pgd,3)
3669 ijump = ni + 2
3670
3671 i = ni
3672 call ngfft( i )
3673 if ( i.ne.ni ) then
3674 write(*,*) 'fft3dvar: NI = ',ni,' I = ',i
3675 call utl_abort('fft3dvar: vector length is not compatible with FFT')
3676 endif
3677 call setfft8(ni)
3678
3679 allocate(pgd2(ni+2,nj))
3680
3681 ! Copy over input data into over-dimensioned array
3682 !$OMP PARALLEL DO PRIVATE (jj,jk,ji,pgd2)
3683 do jk=1,kfield
3684 do jj=1, nj
3685 do ji=1,ni
3686 pgd2(ji,jj)=pgd(ji,jj,jk)
3687 enddo
3688 do ji=ni+1,ni+2
3689 pgd2(ji,jj)=0.0d0
3690 enddo
3691 enddo
3692
3693 call ffft8(pgd2(1,1),1,ijump,nj,kdir)
3694 !* subroutine ffft8( a, inc, jump, lot, isign )
3695 !* a is the array containing input & output data
3696 !* inc is the increment within each data 'vector'
3697 !* (e.g. inc=1 for consecutively stored data)
3698 !* jump is the increment between the start of each data vector
3699 !* lot is the number of data vectors
3700 !* isign = +1 for transform from spectral to gridpoint
3701 !* = -1 for transform from gridpoint to spectral
3702
3703 do jj=1,nj
3704 do ji=1,ni
3705 pgd(ji,jj,jk)=pgd2(ji,jj)
3706 enddo
3707 enddo
3708 enddo
3709 !$OMP END PARALLEL DO
3710
3711 deallocate(pgd2)
3712
3713 end subroutine fft3dvar
3714
3715 subroutine fft3dvar_kij(pgd,kdir)
3716 implicit none
3717
3718 ! Arguments:
3719 real(8) :: pgd(:,:,:)
3720 integer :: kdir
3721
3722 ! Locals:
3723 integer :: kfield,ni,nj
3724 real(8), allocatable :: pgd2(:,:)
3725 integer :: ji,jj,jk,inc,ijump,i
3726
3727 kfield = size(pgd,1)
3728 ni = size(pgd,2)
3729 nj = size(pgd,3)
3730
3731 i = ni
3732 call ngfft( i )
3733 if ( i.ne.ni ) then
3734 write(*,*) 'fft3dvar: NI = ',ni,' I = ',i
3735 call utl_abort('fft3dvar: vector length is not compatible with FFT')
3736 endif
3737 call setfft8(ni)
3738
3739 inc = kfield
3740 ijump = 1
3741 allocate(pgd2(kfield,ni+2))
3742
3743 ! Copy over input data into over-dimensioned array
3744 !$OMP PARALLEL DO PRIVATE (jj,jk,ji,pgd2)
3745 do jj=1, nj
3746 do ji=1,ni
3747 do jk=1,kfield
3748 pgd2(jk,ji)=pgd(jk,ji,jj)
3749 enddo
3750 enddo
3751 do ji=ni+1,ni+2
3752 do jk=1,kfield
3753 pgd2(jk,ji)=0.0d0
3754 enddo
3755 enddo
3756
3757 call ffft8(pgd2(1,1),inc,ijump,kfield,kdir)
3758 !* subroutine ffft8( a, inc, jump, lot, isign )
3759 !* a is the array containing input & output data
3760 !* inc is the increment within each data 'vector'
3761 !* (e.g. inc=1 for consecutively stored data)
3762 !* jump is the increment between the start of each data vector
3763 !* lot is the number of data vectors
3764 !* isign = +1 for transform from spectral to gridpoint
3765 !* = -1 for transform from gridpoint to spectral
3766
3767 do ji=1,ni
3768 do jk=1,kfield
3769 pgd(jk,ji,jj)=pgd2(jk,ji)
3770 enddo
3771 enddo
3772 enddo
3773 !$OMP END PARALLEL DO
3774
3775 deallocate(pgd2)
3776
3777 end subroutine fft3dvar_kij
3778
3779 subroutine ngfft( n )
3780 implicit none
3781
3782 ! Arguments:
3783 integer :: n
3784
3785 ! Locals:
3786 integer, parameter :: l = 3
3787 integer :: k(l), m
3788 data m , k / 8 , 2 , 3 , 5 /
3789 integer :: i,j
3790
3791 if ( n.le.m ) n = m + 1
3792 n = n - 1
37931 n = n + 1
3794 i = n
37952 do 3 j=1,l
3796 if( mod(i,k(j)) .eq. 0 ) go to 4
37973 continue
3798 go to 1
37994 i = i/k(j)
3800 if( i .ne. 1 ) go to 2
3801
3802 end subroutine ngfft
3803
3804
3805end module globalSpectralTransform_mod
