1MODULE advection_mod
2 ! MODULE advection_mod (prefix='adv' category='4. Data Object transformations')
3 !
4 !:Purpose: To perform forward and/or backward advection (based on
5 ! semi-lagrangian trajectories) for both gridStateVector and
6 ! ensemble of gridStateVectors
7 !
8 use ramDisk_mod
9 use midasMpi_mod
10 use mathPhysConstants_mod
11 use earthConstants_mod
12 use ensembleStateVector_mod
13 use gridStateVector_mod
14 use gridStateVectorFileIO_mod
15 use horizontalCoord_mod
16 use verticalCoord_mod
17 use utilities_mod
18 use varNameList_mod
19 implicit none
20 save
21 private
22
23 ! public procedures
24 public :: struct_adv, adv_Setup
25 public :: adv_ensemble_tl, adv_ensemble_ad
26 public :: adv_statevector_tl, adv_statevector_ad
27
28 type :: struct_adv_LevType
29 integer, allocatable :: lonIndex (:,:,:) ! lon, lat, lev
30 integer, allocatable :: latIndex (:,:,:)
31 real(8), allocatable :: interpWeight_BL(:,:,:)
32 real(8), allocatable :: interpWeight_BR(:,:,:)
33 real(8), allocatable :: interpWeight_TL(:,:,:)
34 real(8), allocatable :: interpWeight_TR(:,:,:)
35 end type struct_adv_LevType
36
37 type :: struct_adv_timeStep
38 type(struct_adv_LevType), allocatable :: levType(:)
39 end type struct_adv_timeStep
40
41 type :: struct_adv
42 private
43 integer :: nLev_M
44 integer :: nLev_T
45 integer :: ni, nj
46 integer :: lonPerPE, lonPerPEmax, myLonBeg, myLonEnd
47 integer :: latPerPE, latPerPEmax, myLatBeg, myLatEnd
48 integer, allocatable :: allLonBeg(:), allLatBeg(:)
49 integer :: nTimeStep
50 integer :: timeStepIndexMainSource
51 integer, allocatable :: timeStepIndexSource(:)
52 logical :: singleTimeStepIndexSource
53 type(struct_adv_timeStep), allocatable :: timeStep(:)
54 end type struct_adv
55
56 integer, external :: get_max_rss
57
58 integer, parameter :: MMindex = 1
59 integer, parameter :: THindex = 2
60 integer, parameter :: SFindex = 3
61
62 real(8), parameter :: numGridPts = 1.0d0 ! used to compute numSubStep
63 real(8), parameter :: latitudePatch = 80.0d0 ! this defines latitude where rotated grid used
64
65 integer :: numStepSteeringFlow
66 integer :: lonPerPE, latPerPE
67 integer, allocatable :: allLonBeg(:), allLatBeg(:)
68 integer, allocatable :: numSubStep(:)
69
70 real(8), pointer :: uu_steeringFlow_ptr4d(:,:,:,:)
71 real(8), pointer :: vv_steeringFlow_ptr4d(:,:,:,:)
72 real(8), allocatable :: uu_steeringFlow_mpiGlobal(:,:,:)
73 real(8), allocatable :: vv_steeringFlow_mpiGlobal(:,:,:)
74 real(8), allocatable :: uu_steeringFlow_ThermoLevel(:,:)
75 real(8), allocatable :: vv_steeringFlow_ThermoLevel(:,:)
76
77 real(8), allocatable :: steeringFlowFactor(:)
78
79 real(8) :: steeringFlowDelTsec
80
81 type(struct_hco), pointer :: hco
82 type(struct_vco), pointer :: vco
83
84 logical :: nlat_equalAcrossMpiTasks, nlon_equalAcrossMpiTasks
85
86 ! Control parameter for the level of listing output
87 logical, parameter :: verbose = .false.
88
89CONTAINS
90
91 !--------------------------------------------------------------------------
92 ! adv_Setup
93 !--------------------------------------------------------------------------
94 SUBROUTINE adv_setup(adv, mode, hco_in, vco_in, numStepAdvectedField, &
95 dateStampListAdvectedField, numStepSteeringFlow_in, steeringFlowDelThour, &
96 steeringFlowFactor_in, levTypeList, steeringFlowFilename_opt, &
97 statevector_steeringFlow_opt)
98 implicit none
99
100 ! Arguments:
101 type(struct_adv), intent(inout) :: adv
102 type(struct_hco), pointer, intent(in) :: hco_in
103 type(struct_vco), pointer, intent(in) :: vco_in
104 character(len=*), intent(in) :: mode
105 character(len=*), intent(in) :: levTypeList
106 character(len=*), optional, intent(in) :: steeringFlowFilename_opt
107 integer, intent(in) :: numStepAdvectedField
108 integer, intent(in) :: numStepSteeringFlow_in
109 integer, intent(in) :: dateStampListAdvectedField(numStepAdvectedField)
110 real(8), intent(in) :: steeringFlowFactor_in(vco_in%nLev_M)
111 real(8), intent(in) :: steeringFlowDelThour
112 type(struct_gsv), optional, intent(inout) :: statevector_steeringFlow_opt
113
114 ! Locals:
115 integer :: latIndex0, lonIndex0, latIndex, lonIndex, levIndex, stepIndexSF, stepIndexAF, ierr
116 integer :: levIndexBelow, levIndexAbove
117 integer :: gdxyfll, gdllfxy
118 integer :: nLevType
119 integer, allocatable :: dateStampListSteeringFlow(:)
120 integer, allocatable :: advectedFieldAssociatedStepIndexSF(:)
121 integer, allocatable :: advectionSteeringFlowStartingStepIndex(:)
122 integer, allocatable :: advectionSteeringFlowEndingStepIndex (:)
123 real(8) :: interpWeight_BL, interpWeight_BR, interpWeight_TL, interpWeight_TR
124 real(8), allocatable :: uu_steeringFlow_mpiGlobalTiles(:,:,:,:)
125 real(8), allocatable :: vv_steeringFlow_mpiGlobalTiles(:,:,:,:)
126 real(4) :: xpos_r4, ypos_r4, xposTH_r4, yposTH_r4
127 real(4) :: lonMMbelow_deg_r4, lonMMabove_deg_r4, latMMbelow_deg_r4, latMMabove_deg_r4, lonTH_deg_r4, latTH_deg_r4
128 real(4), allocatable :: xposMM_r4(:,:,:,:), yposMM_r4(:,:,:,:)
129 character(len=1024) :: filename
130
131 type(struct_gsv) :: statevector_steeringFlow
132 logical :: AdvectFileExists
133 integer :: nLev, levTypeIndex, stepIndexSF_start, stepIndexSF_end
134 integer :: myLonBeg, myLonEnd
135 integer :: myLatBeg, myLatEnd
136
137 !
138 !- 1. Set low-level variables
139 !
140 numStepSteeringFlow = numStepSteeringFlow_in
141 adv%nTimeStep = numStepAdvectedField
142
143 allocate(steeringFlowFactor(vco_in%nLev_M))
144 do levIndex = 1, vco_in%nLev_M
145 steeringFlowFactor(levIndex) = steeringFlowFactor_in(levIndex)
146 write(*,*) 'adv_setup: steeringFlowFactor = ', levIndex, steeringFlowFactor(levIndex)
147 end do
148
149 allocate(adv%timeStepIndexSource(numStepAdvectedField))
150
151 if (vco_in%Vcode /= 5002 .and. vco_in%Vcode /= 5005 ) then
152 call utl_abort('adv_setup: Only vCode 5002 and 5005 are currently supported!')
153 end if
154
155 select case(trim(levTypeList))
156 case ('MMLevsOnly')
157 nLevType = 1
158 case ('allLevs')
159 nLevType = 3
160 case default
161 write(*,*)
162 write(*,*) 'Unsupported levTypeList: ', trim(levTypeList)
163 call utl_abort('adv_setup')
164 end select
165
166 !- 1.1 Mode
167 select case(trim(mode))
168 case ('fromFirstTimeIndex')
169 adv%timeStepIndexMainSource = 1
170 adv%timeStepIndexSource(:) = adv%timeStepIndexMainSource
171 adv%singleTimeStepIndexSource= .true.
172 case ('fromMiddleTimeIndex')
173 if (mod(numStepAdvectedField,2) == 0) then
174 call utl_abort('adv_setup: numStepAdvectedField cannot be even with direction=fromMiddleTimeIndex')
175 end if
176 adv%timeStepIndexMainSource = (numStepAdvectedField+1)/2
177 adv%timeStepIndexSource(:) = adv%timeStepIndexMainSource
178 adv%singleTimeStepIndexSource= .true.
179 case ('fromLastTimeIndex')
180 adv%timeStepIndexMainSource = numStepAdvectedField
181 adv%timeStepIndexSource(:) = adv%timeStepIndexMainSource
182 adv%singleTimeStepIndexSource= .true.
183 case ('towardFirstTimeIndex','towardFirstTimeIndexInverse')
184 adv%timeStepIndexMainSource = 1
185 do stepIndexAF = 1, numStepAdvectedField
186 adv%timeStepIndexSource(stepIndexAF) = stepIndexAF
187 end do
188 adv%singleTimeStepIndexSource = .false.
189 case('towardMiddleTimeIndex','towardMiddleTimeIndexInverse')
190 if (mod(numStepAdvectedField,2) == 0) then
191 call utl_abort('adv_setup: numStepAdvectedField cannot be even with direction=towardMiddleTimeIndex')
192 end if
193 adv%timeStepIndexMainSource = (numStepAdvectedField+1)/2
194 do stepIndexAF = 1, numStepAdvectedField
195 adv%timeStepIndexSource(stepIndexAF) = stepIndexAF
196 end do
197 adv%singleTimeStepIndexSource = .false.
198 case default
199 write(*,*)
200 write(*,*) 'Unsupported mode : ', trim(mode)
201 call utl_abort('adv_setup')
202 end select
203
204 !- Set some important values
205 steeringFlowDelTsec = steeringFlowDelThour*3600.0D0
206
207 !- 1.2 Grid Size
208 hco => hco_in
209 adv%ni = hco%ni
210 adv%nj = hco%nj
211 vco => vco_in
212 adv%nLev_M = vco%nLev_M
213 adv%nLev_T = vco%nLev_T
214
215 call mmpi_setup_latbands(adv%nj, adv%latPerPE, adv%latPerPEmax, adv%myLatBeg, adv%myLatEnd, &
216 divisible_opt=nlat_equalAcrossMpiTasks)
217 call mmpi_setup_lonbands(adv%ni, adv%lonPerPE, adv%lonPerPEmax, adv%myLonBeg, adv%myLonEnd, &
218 divisible_opt=nlon_equalAcrossMpiTasks)
219 allocate(adv%allLonBeg(mmpi_npex))
220 call rpn_comm_allgather(adv%myLonBeg,1,"mpi_integer", &
221 adv%allLonBeg,1,"mpi_integer","EW",ierr)
222 allocate(adv%allLatBeg(mmpi_npey))
223 call rpn_comm_allgather(adv%myLatBeg,1,"mpi_integer", &
224 adv%allLatBeg,1,"mpi_integer","NS",ierr)
225
226 lonPerPE = adv%lonPerPE
227 latPerPE = adv%latPerPE
228 allocate(allLonBeg(mmpi_npex))
229 allocate(allLatBeg(mmpi_npey))
230 allLonBeg(:) = adv%allLonBeg(:)
231 allLatBeg(:) = adv%allLatBeg(:)
232
233 !- 1.3 Memory allocation
234 myLonBeg = adv%myLonBeg
235 myLonEnd = adv%myLonEnd
236 myLatBeg = adv%myLatBeg
237 myLatEnd = adv%myLatEnd
238
239 nLev = adv%nLev_M
240
241 allocate(adv%timeStep(numStepAdvectedField))
242
243 do stepIndexAF = 1, numStepAdvectedField
244
245 if (stepIndexAF == adv%timeStepIndexMainSource) cycle ! no interpolation needed for this time step
246
247 allocate(adv%timeStep(stepIndexAF)%levType(nLevType))
248 do levTypeIndex = 1, nLevType ! ( 1=MM, 2=TH, 3=SF )
249 if (levTypeIndex == MMindex ) then
250 nLev = adv%nLev_M
251 else if (levTypeIndex == THindex) then
252 nLev = adv%nLev_T
253 else
254 nLev = 1
255 end if
256 allocate(adv%timeStep(stepIndexAF)%levType(levTypeIndex)%lonIndex (myLonBeg:myLonEnd,myLatBeg:myLatEnd,nLev))
257 allocate(adv%timeStep(stepIndexAF)%levType(levTypeIndex)%latIndex (myLonBeg:myLonEnd,myLatBeg:myLatEnd,nLev))
258 allocate(adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BL(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nLev))
259 allocate(adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BR(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nLev))
260 allocate(adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TL(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nLev))
261 allocate(adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TR(myLonBeg:myLonEnd,myLatBeg:myLatEnd,nLev))
262 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BL(:,:,:) = 1.0d0
263 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BR(:,:,:) = 0.0d0
264 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TL(:,:,:) = 0.0d0
265 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TR(:,:,:) = 0.0d0
266 end do
267 end do
268
269 !
270 !- 2. Read in the wind data to use for advection
271 !
272 allocate(dateStampListSteeringFlow(numStepSteeringFlow))
273
274 if (present(steeringFlowFilename_opt) ) then
275
276 if (mmpi_myid == 0) then
277 write(*,*)
278 write(*,*) 'steeringFlow source taken from input file = ', trim(steeringFlowFilename_opt)
279 end if
280
281 do stepIndexSF = 1, numStepSteeringFlow
282 call incdatr(dateStampListSteeringFlow(stepIndexSF), dateStampListAdvectedField(1), &
283 real(stepIndexSF-1,8)*steeringFlowDelThour)
284 end do
285
286 call gsv_allocate(statevector_steeringFlow,numStepSteeringFlow, hco, vco, &
287 dateStampList_opt=dateStampListSteeringFlow, &
288 varNames_opt=(/'UU','VV','P0'/), mpi_local_opt=.true., &
289 hInterpolateDegree_opt='LINEAR')
290
291 fileName = ram_fullWorkingPath(trim(steeringFlowFilename_opt))
292 inquire(file=trim(fileName),exist=AdvectFileExists)
293 write(*,*) 'AdvectFileExists', AdvectFileExists
294 do stepIndexSF = 1, numStepSteeringFlow
295 call gio_readFromFile( statevector_steeringFlow, fileName, ' ', ' ', stepIndex_opt=stepIndexSF, &
296 containsFullField_opt=.true.)
297 end do
298
299 ierr = ram_remove(fileName)
300
301 call gsv_getField(statevector_steeringFlow, uu_steeringFlow_ptr4d, 'UU')
302 call gsv_getField(statevector_steeringFlow, vv_steeringFlow_ptr4d, 'VV')
303
304 else if (present(statevector_steeringFlow_opt)) then
305
306 if (mmpi_myid == 0) then
307 write(*,*)
308 write(*,*) 'steeringFlow source = input gridStateVector'
309 write(*,*) numStepSteeringFlow
310 write(*,*) statevector_steeringFlow_opt%dateStampList(:)
311 end if
312 dateStampListSteeringFlow(:) = statevector_steeringFlow_opt%dateStampList(:)
313 call gsv_getField(statevector_steeringFlow_opt, uu_steeringFlow_ptr4d, 'UU')
314 call gsv_getField(statevector_steeringFlow_opt, vv_steeringFlow_ptr4d, 'VV')
315
316 else
317 call utl_abort('adv_setup: steeringFlow source was not provided!')
318 end if
319
320 !
321 !- 3. Advection setup
322 !
323
324 !- 3.1 Match the stepIndex between the reference flow and the fields to be advected
325 allocate(advectedFieldAssociatedStepIndexSF(numStepAdvectedField))
326 advectedFieldAssociatedStepIndexSF(:) = -1
327 do stepIndexAF = 1, numStepAdvectedField
328 do stepIndexSF = 1, numStepSteeringFlow
329 if ( dateStampListAdvectedField(stepIndexAF) == dateStampListSteeringFlow(stepIndexSF) ) then
330 advectedFieldAssociatedStepIndexSF(stepIndexAF) = stepIndexSF
331 if (mmpi_myid == 0) then
332 write(*,*)
333 write(*,*) 'stepIndex Match', stepIndexAF, stepIndexSF
334 end if
335 exit
336 end if
337 end do
338 if ( advectedFieldAssociatedStepIndexSF(stepIndexAF) == -1 ) then
339 call utl_abort('adv_setup: no match between dateStampListAdvectedField and dateStampListSteeringFlow')
340 end if
341 end do
342
343 !- 3.2 Set starting, ending and direction parameters
344 allocate(advectionSteeringFlowStartingStepIndex(numStepAdvectedField))
345 allocate(advectionSteeringFlowEndingStepIndex (numStepAdvectedField))
346
347 select case(trim(mode))
348 case ('fromFirstTimeIndex','fromMiddleTimeIndex','fromLastTimeIndex', &
349 'towardFirstTimeIndexInverse','towardMiddleTimeIndexInverse')
350 do stepIndexAF = 1, numStepAdvectedField
351 advectionSteeringFlowStartingStepIndex(stepIndexAF) = advectedFieldAssociatedStepIndexSF(adv%timeStepIndexMainSource)
352 advectionSteeringFlowEndingStepIndex (stepIndexAF) = advectedFieldAssociatedStepIndexSF(stepIndexAF)
353 end do
354 case ('towardFirstTimeIndex','towardMiddleTimeIndex')
355 do stepIndexAF = 1, numStepAdvectedField
356 advectionSteeringFlowStartingStepIndex(stepIndexAF) = advectedFieldAssociatedStepIndexSF(stepIndexAF)
357 advectionSteeringFlowEndingStepIndex (stepIndexAF) = advectedFieldAssociatedStepIndexSF(adv%timeStepIndexMainSource)
358 end do
359 case default
360 write(*,*)
361 write(*,*) 'Oops! This should never happen. Check the code...'
362 call utl_abort('adv_setup')
363 end select
364
365 !
366 !- 4. Perform the advection (backward and/or forward)
367 !
368 if (mmpi_myid == 0) write(*,*) 'setupAdvectAmplitude: starting'
369 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
370
371 allocate(numSubStep(adv%nj))
372 allocate(uu_steeringFlow_mpiGlobal(numStepSteeringFlow, adv%ni, adv%nj))
373 allocate(vv_steeringFlow_mpiGlobal(numStepSteeringFlow, adv%ni, adv%nj))
374
375 allocate(uu_steeringFlow_mpiGlobalTiles(numStepSteeringFlow, adv%lonPerPE, adv%latPerPE, mmpi_nprocs))
376 allocate(vv_steeringFlow_mpiGlobalTiles(numStepSteeringFlow, adv%lonPerPE, adv%latPerPE, mmpi_nprocs))
377 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
378
379 if ( nLevType >= THindex ) then
380 allocate(xposMM_r4(numStepSteeringFlow,myLonBeg:myLonEnd,myLatBeg:myLatEnd,adv%nLev_M))
381 allocate(yposMM_r4(numStepSteeringFlow,myLonBeg:myLonEnd,myLatBeg:myLatEnd,adv%nLev_M))
382 end if
383
384 !- 4.1 Compute the trajectories on the momentum levels
385 levTypeIndex = MMindex
386 nLev = adv%nLev_M
387
388 do levIndex = 1, nLev ! loop over levels
389
390 if (mmpi_myid == 0) write(*,*) 'setupAdvectAmplitude: levIndex = ', levIndex
391
392 call processSteeringFlow(levTypeIndex, levIndex, & ! IN
393 uu_steeringFlow_mpiGlobalTiles, vv_steeringFlow_mpiGlobalTiles, & ! OUT
394 adv%nLev_M, adv%nLev_T, myLatBeg, myLatEnd) ! IN
395
396 do stepIndexAF = 1, numStepAdvectedField
397 if (stepIndexAF == adv%timeStepIndexMainSource) cycle ! no interpolation needed for this time step
398
399 stepIndexSF_start = advectionSteeringFlowStartingStepIndex(stepIndexAF)
400 stepIndexSF_end = advectionSteeringFlowEndingStepIndex (stepIndexAF)
401
402 ! loop over all initial grid points within tile for determining trajectories
403 do latIndex0 = adv%myLatBeg, adv%myLatEnd
404 do lonIndex0 = adv%myLonBeg, adv%myLonEnd
405
406 call calcTrajectory(xpos_r4, ypos_r4, & ! OUT
407 latIndex0, lonIndex0, levIndex, stepIndexSF_start, stepIndexSF_end) ! IN
408
409 if ( nLevType >= THindex ) then
410 xposMM_r4(stepIndexAF,lonIndex0,latIndex0,levIndex) = xpos_r4
411 yposMM_r4(stepIndexAF,lonIndex0,latIndex0,levIndex) = ypos_r4
412 end if
413
414 call calcWeights(lonIndex, latIndex, interpWeight_BL, interpWeight_BR, & ! OUT
415 interpWeight_TL, interpWeight_TR, & ! OUT
416 xpos_r4, ypos_r4) ! IN
417
418 ! store the final position of the trajectory and interp weights
419 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%lonIndex (lonIndex0,latIndex0,levIndex) = lonIndex
420 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%latIndex (lonIndex0,latIndex0,levIndex) = latIndex
421 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BL(lonIndex0,latIndex0,levIndex) = interpWeight_BL
422 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BR(lonIndex0,latIndex0,levIndex) = interpWeight_BR
423 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TL(lonIndex0,latIndex0,levIndex) = interpWeight_TL
424 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TR(lonIndex0,latIndex0,levIndex) = interpWeight_TR
425
426 end do ! lonIndex0
427 end do ! latIndex0
428
429 end do ! stepIndexAF
430
431 end do ! levIndex
432
433 !- 4.2 Thermodynamic levels: Interpolate vertically the positions found on the momentum levels
434 if ( nLevType >= THindex ) then
435 levTypeIndex = THindex
436 nLev = adv%nLev_T
437
438 do levIndex = 1, nLev ! loop over levels
439 do stepIndexAF = 1, numStepAdvectedField
440 if (stepIndexAF == adv%timeStepIndexMainSource) cycle ! no interpolation needed for this time step
441
442 do latIndex0 = adv%myLatBeg, adv%myLatEnd
443 do lonIndex0 = adv%myLonBeg, adv%myLonEnd
444
445 if (levIndex == 1 .and. vco%Vcode == 5002) then
446 ! use top momentum level amplitudes for top thermo level
447 xposTH_r4 = xposMM_r4(stepIndexAF,lonIndex0,latIndex0,1)
448 yposTH_r4 = yposMM_r4(stepIndexAF,lonIndex0,latIndex0,1)
449 else if (levIndex == nLev) then
450 ! use surface momentum level amplitudes for surface thermo level
451 xposTH_r4 = xposMM_r4(stepIndexAF,lonIndex0,latIndex0,adv%nLev_M)
452 yposTH_r4 = yposMM_r4(stepIndexAF,lonIndex0,latIndex0,adv%nLev_M)
453 else
454 ! for other levels, interpolate momentum positions to get thermo positions (as in GEM)
455 if (vco%Vcode == 5002) then
456 levIndexBelow = levIndex
457 levIndexAbove = levIndex-1
458 else
459 levIndexBelow = levIndex+1
460 levIndexAbove = levIndex
461 end if
462
463 ierr = gdllfxy(hco%EZscintID, latMMbelow_deg_r4, lonMMbelow_deg_r4, &
464 xposMM_r4(stepIndexAF,lonIndex0,latIndex0,levIndexBelow), &
465 yposMM_r4(stepIndexAF,lonIndex0,latIndex0,levIndexBelow), 1)
466
467 ierr = gdllfxy(hco%EZscintID, latMMabove_deg_r4, lonMMabove_deg_r4, &
468 xposMM_r4(stepIndexAF,lonIndex0,latIndex0,levIndexAbove), &
469 yposMM_r4(stepIndexAF,lonIndex0,latIndex0,levIndexAbove), 1)
470
471 if (lonMMbelow_deg_r4 < 0.0) lonMMbelow_deg_r4 = lonMMbelow_deg_r4 + 360.0
472 if (lonMMabove_deg_r4 < 0.0) lonMMabove_deg_r4 = lonMMabove_deg_r4 + 360.0
473
474 if ( abs(lonMMbelow_deg_r4 - lonMMabove_deg_r4) > 180.0 ) then
475 if (lonMMbelow_deg_r4 > 180.0 ) then
476 lonMMbelow_deg_r4 = lonMMbelow_deg_r4 - 360.0
477 else
478 lonMMabove_deg_r4 = lonMMabove_deg_r4 - 360.0
479 end if
480 end if
481 lonTH_deg_r4 = 0.5 * (lonMMbelow_deg_r4 + lonMMabove_deg_r4)
482 if (lonTH_deg_r4 < 0.0) lonTH_deg_r4 = lonTH_deg_r4 + 360.0
483 latTH_deg_r4 = 0.5 * (latMMbelow_deg_r4 + latMMabove_deg_r4)
484
485 ierr = gdxyfll(hco%EZscintID, xposTH_r4, yposTH_r4, &
486 latTH_deg_r4, lonTH_deg_r4, 1)
487 end if
488
489 ! Compute weights
490 call calcWeights(lonIndex, latIndex, interpWeight_BL, interpWeight_BR, & ! OUT
491 interpWeight_TL, interpWeight_TR, & ! OUT
492 xposTH_r4, yposTH_r4) ! IN
493
494 ! Store the final position of the trajectory and interp weights
495 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%lonIndex (lonIndex0,latIndex0,levIndex) = lonIndex
496 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%latIndex (lonIndex0,latIndex0,levIndex) = latIndex
497 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BL(lonIndex0,latIndex0,levIndex) = interpWeight_BL
498 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BR(lonIndex0,latIndex0,levIndex) = interpWeight_BR
499 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TL(lonIndex0,latIndex0,levIndex) = interpWeight_TL
500 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TR(lonIndex0,latIndex0,levIndex) = interpWeight_TR
501 end do
502 end do
503 end do
504 end do
505
506 deallocate(xposMM_r4)
507 deallocate(yposMM_r4)
508
509 end if
510
511 !- 4.3 Surface level: Use the positions from the lowest momentum levels
512 if ( nLevType == SFindex ) then
513 do stepIndexAF = 1, numStepAdvectedField
514 if (stepIndexAF == adv%timeStepIndexMainSource) cycle ! no interpolation needed for this time step
515 adv%timeStep(stepIndexAF)%levType(SFindex)%lonIndex (:,:,1) = &
516 adv%timeStep(stepIndexAF)%levType(MMindex)%lonIndex (:,:,adv%nLev_M)
517 adv%timeStep(stepIndexAF)%levType(SFindex)%latIndex (:,:,1) = &
518 adv%timeStep(stepIndexAF)%levType(MMindex)%latIndex (:,:,adv%nLev_M)
519 adv%timeStep(stepIndexAF)%levType(SFindex)%interpWeight_BL(:,:,1) = &
520 adv%timeStep(stepIndexAF)%levType(MMindex)%interpWeight_BL(:,:,adv%nLev_M)
521 adv%timeStep(stepIndexAF)%levType(SFindex)%interpWeight_BR(:,:,1) = &
522 adv%timeStep(stepIndexAF)%levType(MMindex)%interpWeight_BR(:,:,adv%nLev_M)
523 adv%timeStep(stepIndexAF)%levType(SFindex)%interpWeight_TL(:,:,1) = &
524 adv%timeStep(stepIndexAF)%levType(MMindex)%interpWeight_TL(:,:,adv%nLev_M)
525 adv%timeStep(stepIndexAF)%levType(SFindex)%interpWeight_TR(:,:,1) = &
526 adv%timeStep(stepIndexAF)%levType(MMindex)%interpWeight_TR(:,:,adv%nLev_M)
527 end do
528 end if
529
530 deallocate(numSubStep)
531 deallocate(allLonBeg)
532 deallocate(allLatBeg)
533 deallocate(uu_steeringFlow_mpiGlobalTiles)
534 deallocate(vv_steeringFlow_mpiGlobalTiles)
535 deallocate(uu_steeringFlow_mpiGlobal)
536 deallocate(vv_steeringFlow_mpiGlobal)
537 if (present(steeringFlowFilename_opt) ) then
538 call gsv_deallocate(statevector_steeringFlow)
539 end if
540
541 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
542 if (mmpi_myid == 0) write(*,*) 'adv_setup: done'
543
544 end SUBROUTINE adv_setup
545
546 !--------------------------------------------------------------------------
547 ! processSteeringFlow
548 !--------------------------------------------------------------------------
549 SUBROUTINE processSteeringFlow (levTypeIndex, levIndex, &
550 uu_steeringFlow_mpiGlobalTiles, vv_steeringFlow_mpiGlobalTiles, &
551 nLev_M, nLev_T, myLatBeg, myLatEnd)
552 implicit none
553
554 ! Arguments:
555 integer, intent(in) :: levTypeIndex
556 integer, intent(in) :: levIndex
557 integer, intent(in) :: nLev_M
558 integer, intent(in) :: nLev_T
559 integer, intent(in) :: myLatBeg
560 integer, intent(in) :: myLatEnd
561 real(8), intent(out) :: uu_steeringFlow_mpiGlobalTiles(:,:,:,:)
562 real(8), intent(out) :: vv_steeringFlow_mpiGlobalTiles(:,:,:,:)
563
564 ! Locals:
565 integer :: stepIndexSF, nsize, ierr
566 integer :: procID, procIDx, procIDy, lonIndex, latIndex
567 integer :: lonIndex_mpiglobal, latIndex_mpiglobal
568 integer :: levIndexBelow, levIndexAbove
569 real(8) :: uu, vv, latAdvect
570
571 nsize = lonPerPE*latPerPE
572
573 if ( levTypeIndex == MMindex ) then
574 ! No vertical interpolation is needed
575 do stepIndexSF = 1, numStepSteeringFlow
576 ! gather the winds for this level
577 call rpn_comm_allgather(uu_steeringFlow_ptr4d(:,:,levIndex,stepIndexSF) , nsize, "mpi_double_precision", &
578 uu_steeringFlow_mpiGlobalTiles(stepIndexSF,:,:,:), nsize, "mpi_double_precision", &
579 "GRID", ierr )
580 call rpn_comm_allgather(vv_steeringFlow_ptr4d(:,:,levIndex,stepIndexSF) , nsize, "mpi_double_precision", &
581 vv_steeringFlow_mpiGlobalTiles(stepIndexSF,:,:,:), nsize, "mpi_double_precision", &
582 "GRID", ierr )
583 end do
584 else if (levTypeIndex == THindex ) then
585 ! Vertical interpolation is needed...
586 ! The adopted approach follows the vertical interpolation for the amplitude fields in bMatrixEnsemble_mod
587 do stepIndexSF = 1, numStepSteeringFlow
588 if (levIndex == 1 .and. vco%Vcode == 5002) then
589 ! use top momentum level amplitudes for top thermo level
590 uu_steeringFlow_ThermoLevel(:,:) = uu_steeringFlow_ptr4d(:,:,1,stepIndexSF)
591 vv_steeringFlow_ThermoLevel(:,:) = vv_steeringFlow_ptr4d(:,:,1,stepIndexSF)
592 else if (levIndex == nLev_T) then
593 ! use surface momentum level amplitudes for surface thermo level
594 uu_steeringFlow_ThermoLevel(:,:) = uu_steeringFlow_ptr4d(:,:,nLev_M,stepIndexSF)
595 vv_steeringFlow_ThermoLevel(:,:) = vv_steeringFlow_ptr4d(:,:,nLev_M,stepIndexSF)
596 else
597 ! for other levels, interpolate momentum winds to get thermo winds
598 if (vco%Vcode == 5002) then
599 levIndexBelow = levIndex
600 levIndexAbove = levIndex-1
601 else
602 levIndexBelow = levIndex+1
603 levIndexAbove = levIndex
604 end if
605 !$OMP PARALLEL DO PRIVATE (latIndex)
606 do latIndex = myLatBeg, myLatEnd
607 uu_steeringFlow_ThermoLevel(:,latIndex) = 0.5d0*( uu_steeringFlow_ptr4d(:,latIndex,levIndexAbove,stepIndexSF) + &
608 uu_steeringFlow_ptr4d(:,latIndex,levIndexBelow,stepIndexSF) )
609 vv_steeringFlow_ThermoLevel(:,latIndex) = 0.5d0*( vv_steeringFlow_ptr4d(:,latIndex,levIndexAbove,stepIndexSF) + &
610 vv_steeringFlow_ptr4d(:,latIndex,levIndexBelow,stepIndexSF) )
611 end do
612 !$OMP END PARALLEL DO
613 end if
614
615 ! gather the INTERPOLATED winds for this level
616 call rpn_comm_allgather(uu_steeringFlow_ThermoLevel , nsize, "mpi_double_precision", &
617 uu_steeringFlow_mpiGlobalTiles(stepIndexSF,:,:,:), nsize, "mpi_double_precision", &
618 "GRID", ierr )
619 call rpn_comm_allgather(uu_steeringFlow_ThermoLevel , nsize, "mpi_double_precision", &
620 vv_steeringFlow_mpiGlobalTiles(stepIndexSF,:,:,:), nsize, "mpi_double_precision", &
621 "GRID", ierr )
622
623 end do
624
625 else
626 call utl_abort('processSteeringFlow: invalid levTypeIndex')
627 end if
628
629 ! rearrange gathered winds for convenience
630 do procIDy = 0, (mmpi_npey-1)
631 do procIDx = 0, (mmpi_npex-1)
632 procID = procIDx + procIDy*mmpi_npex
633 do latIndex = 1, latPerPE
634 latIndex_mpiglobal = latIndex + allLatBeg(procIDy+1) - 1
635 do lonIndex = 1, lonPerPE
636 lonIndex_mpiglobal = lonIndex + allLonBeg(procIDx+1) - 1
637 uu_steeringFlow_mpiGlobal(:, lonIndex_mpiglobal, latIndex_mpiglobal) = uu_steeringFlow_mpiGlobalTiles(:, lonIndex, latIndex, procID+1)
638 vv_steeringFlow_mpiGlobal(:, lonIndex_mpiglobal, latIndex_mpiglobal) = vv_steeringFlow_mpiGlobalTiles(:, lonIndex, latIndex, procID+1)
639 end do
640 end do
641 end do
642 end do
643
644 ! determine the number of time steps required as a function of latitude
645 do latIndex = 1, hco%nj
646 latAdvect = hco%lat(latIndex)
647 if (abs(latAdvect) < latitudePatch*MPC_RADIANS_PER_DEGREE_R8) then
648 uu = maxval(abs(uu_steeringFlow_mpiGlobal(:,:,latIndex) /(ec_ra*cos(latAdvect)))) ! in rad/s
649 vv = maxval(abs(vv_steeringFlow_mpiGlobal(:,:,latIndex) / ec_ra)) ! in rad/s
650 else
651 uu = maxval(abs(uu_steeringFlow_mpiGlobal(:,:,latIndex) / ec_ra)) ! in rad/s
652 vv = maxval(abs(vv_steeringFlow_mpiGlobal(:,:,latIndex) / ec_ra)) ! in rad/s
653 end if
654 numSubStep(latIndex) = max( 1, &
655 nint( (steeringFlowDelTsec * steeringFlowFactor(levIndex) * uu) / (numGridPts*(hco%lon(2)-hco%lon(1))) ), &
656 nint( (steeringFlowDelTsec * steeringFlowFactor(levIndex) * vv) / (numGridPts*(hco%lat(2)-hco%lat(1))) ) )
657 end do
658 if (mmpi_myid == 0) write(*,*) 'min and max of numSubStep',minval(numSubStep(:)),maxval(numSubStep(:))
659
660 end SUBROUTINE processSteeringFlow
661
662 !--------------------------------------------------------------------------
663 ! calcTrajectory
664 !--------------------------------------------------------------------------
665 SUBROUTINE calcTrajectory(xpos_r4, ypos_r4, latIndex0, lonIndex0, &
666 levIndex, stepIndexSF_start, stepIndexSF_end)
667 implicit none
668
669 ! Arguments:
670 real(4), intent(out) :: xpos_r4
671 real(4), intent(out) :: ypos_r4
672 integer, intent(in) :: latIndex0
673 integer, intent(in) :: lonIndex0
674 integer, intent(in) :: stepIndexSF_start
675 integer, intent(in) :: stepIndexSF_end
676 integer, intent(in) :: levIndex
677
678 ! Locals:
679 integer :: subStepIndex, stepIndexSF, ierr, gdxyfll, latIndex, lonIndex
680 integer :: alfa, ni, nj, stepIndex_direction, stepIndex_first, stepIndex_last
681 real(8) :: uu, vv, subDelT, lonAdvect, latAdvect
682 real(8) :: uu_p, vv_p, lonAdvect_p, latAdvect_p, Gcoef, Scoef
683 real(4) :: lonAdvect_deg_r4, latAdvect_deg_r4
684
685 ni = hco%ni
686 nj = hco%nj
687
688 subDelT = steeringFlowDelTsec/real(numSubStep(latIndex0),8) ! in seconds
689
690 ! position at the initial time of back trajectory
691 lonAdvect = hco%lon(lonIndex0) ! in radians
692 latAdvect = hco%lat(latIndex0)
693 lonIndex = lonIndex0 ! index
694 latIndex = latIndex0
695 xpos_r4 = real(lonIndex,4)
696 ypos_r4 = real(latIndex,4)
697
698 ! initial positions in rotated coordinate system
699 lonAdvect_p = 0.0d0
700 latAdvect_p = 0.0d0
701
702 if (verbose) then
703 write(*,*) 'final lonAdvect,latAdvect=', &
704 lonAdvect*MPC_DEGREES_PER_RADIAN_R8, &
705 latAdvect*MPC_DEGREES_PER_RADIAN_R8
706 write(*,*) 'numSubStep=', numSubStep(latIndex0)
707 end if
708
709 ! time stepping strategy
710 if (stepIndexSF_end > stepIndexSF_start) then
711 ! back trajectory , stepping backwards
712 stepIndex_first = stepIndexSF_end-1
713 stepIndex_last = stepIndexSF_start
714 stepIndex_direction = -1
715 else if (stepIndexSF_end < stepIndexSF_start) then
716 ! forward trajectory, stepping forward
717 stepIndex_first = stepIndexSF_end
718 stepIndex_last = stepIndexSF_start-1
719 stepIndex_direction = 1
720 else
721 call utl_abort('calcTrajAndWeights: fatal error with stepIndexSF')
722 end if
723
724 do stepIndexSF = stepIndex_first, stepIndex_last, stepIndex_direction
725
726 if (verbose) write(*,*) 'stepIndexSF,lonIndex,latIndex=',stepIndexSF,lonIndex,latIndex
727
728 do subStepIndex = 1, numSubStep(latIndex0)
729
730 alfa = (subStepIndex-1)/numSubStep(latIndex0)
731 ! perform one timestep
732 if (abs(hco%lat(latIndex0)) < latitudePatch*MPC_RADIANS_PER_DEGREE_R8) then
733 ! points away from pole, handled normally
734 ! determine wind at current location (now at BL point)
735 uu = ( alfa *uu_steeringFlow_mpiGlobal(stepIndexSF ,lonIndex,latIndex) + &
736 (1-alfa)*uu_steeringFlow_mpiGlobal(stepIndexSF+1,lonIndex,latIndex) ) &
737 /(ec_ra*cos(hco%lat(latIndex))) ! in rad/s
738 vv = ( alfa*vv_steeringFlow_mpiGlobal(stepIndexSF ,lonIndex,latIndex) + &
739 (1-alfa)*vv_steeringFlow_mpiGlobal(stepIndexSF+1,lonIndex,latIndex) ) &
740 /ec_ra
741 ! apply user-specified scale factor to advecting winds
742 uu = steeringFlowFactor(levIndex) * uu
743 vv = steeringFlowFactor(levIndex) * vv
744
745 ! compute next position
746 lonAdvect = lonAdvect + real(stepIndex_direction,8)*subDelT*uu ! in radians
747 latAdvect = latAdvect + real(stepIndex_direction,8)*subDelT*vv
748
749 if (verbose) then
750 write(*,*) 'not near pole, lonAdvect,latAdvect,uu,vv=', &
751 lonAdvect*MPC_DEGREES_PER_RADIAN_R8, &
752 latAdvect*MPC_DEGREES_PER_RADIAN_R8,uu,vv
753 end if
754 else
755 ! points near pole, handled in a special way
756 ! determine wind at current location (now at BL point)
757 uu = alfa *uu_steeringFlow_mpiGlobal(stepIndexSF ,lonIndex,latIndex) + &
758 (1-alfa)*uu_steeringFlow_mpiGlobal(stepIndexSF+1,lonIndex,latIndex) ! in m/s
759 vv = alfa *vv_steeringFlow_mpiGlobal(stepIndexSF ,lonIndex,latIndex) + &
760 (1-alfa)*vv_steeringFlow_mpiGlobal(stepIndexSF+1,lonIndex,latIndex)
761 ! transform wind vector into rotated coordinate system
762 Gcoef = ( cos(latAdvect)*cos(hco%lat(latIndex0)) + &
763 sin(latAdvect)*sin(hco%lat(latIndex0))*cos(lonAdvect-hco%lon(lonIndex0)) ) / &
764 cos(latAdvect_p)
765 Scoef = ( sin(hco%lat(latIndex0))*sin(lonAdvect-hco%lon(lonIndex0)) ) / &
766 cos(latAdvect_p)
767 uu_p = Gcoef * uu - Scoef * vv ! in m/s
768 vv_p = Scoef * uu + Gcoef * vv
769
770 ! apply user-specified scale factor to advecting winds
771 uu_p = steeringFlowFactor(levIndex) * uu_p ! in m/s
772 vv_p = steeringFlowFactor(levIndex) * vv_p
773
774 ! compute next position (in rotated coord system)
775 lonAdvect_p = lonAdvect_p + real(stepIndex_direction,8)*subDelT*uu_p/(ec_ra*cos(latAdvect_p)) ! in radians
776 latAdvect_p = latAdvect_p + real(stepIndex_direction,8)*subDelT*vv_p/ec_ra
777
778 if (verbose) then
779 write(*,*) ' near pole, uu_p,vv_p,Gcoef,Scoef=', &
780 uu_p, vv_p, Gcoef, Scoef
781 write(*,*) ' near pole, lonAdvect_p,latAdvect_p=', &
782 lonAdvect_p*MPC_DEGREES_PER_RADIAN_R8, &
783 latAdvect_p*MPC_DEGREES_PER_RADIAN_R8
784 end if
785
786 ! compute lon/lat in original coordinate system
787 lonAdvect = hco%lon(lonIndex0) + &
788 atan2( cos(latAdvect_p)*sin(lonAdvect_p) , &
789 ( cos(latAdvect_p)*cos(lonAdvect_p)*cos(hco%lat(latIndex0)) - &
790 sin(latAdvect_p)*sin(hco%lat(latIndex0)) ) )
791 latAdvect = asin( cos(latAdvect_p)*cos(lonAdvect_p)*sin(hco%lat(latIndex0)) + &
792 sin(latAdvect_p)*cos(hco%lat(latIndex0)) )
793
794 if (verbose) then
795 write(*,*) ' near pole, lonAdvect,latAdvect=', &
796 lonAdvect*MPC_DEGREES_PER_RADIAN_R8, &
797 latAdvect*MPC_DEGREES_PER_RADIAN_R8
798 end if
799 end if
800
801 ! convert lon/lat position into index
802 lonAdvect_deg_r4 = real(lonAdvect,4)* MPC_DEGREES_PER_RADIAN_R4
803 latAdvect_deg_r4 = real(latAdvect,4)* MPC_DEGREES_PER_RADIAN_R4
804 ierr = gdxyfll(hco%EZscintID, xpos_r4, ypos_r4, &
805 latAdvect_deg_r4, lonAdvect_deg_r4, 1)
806
807 ! determine the bottom-left grid point
808 lonIndex = floor(xpos_r4)
809 latIndex = floor(ypos_r4)
810
811 ! check if position is east of the grid
812 if (floor(xpos_r4) > ni) then
813 if (verbose) then
814 write(*,*) 'lonIndex0,latIndex0,lon,lat,stepIndexSF,x,y xpos_r4 > ni :', &
815 lonIndex0,latIndex0,lonAdvect*MPC_DEGREES_PER_RADIAN_R8, &
816 latAdvect*MPC_DEGREES_PER_RADIAN_R8,stepIndexSF,xpos_r4,ypos_r4
817 end if
818 ! add 10*epsilon(real*4) to ensure do not go too far due to limited precision
819 lonAdvect = lonAdvect - 2.0D0*MPC_PI_R8 + 10.0D0*real(epsilon(1.0),8)
820 lonAdvect_deg_r4 = real(lonAdvect,4)* MPC_DEGREES_PER_RADIAN_R4
821 latAdvect_deg_r4 = real(latAdvect,4)* MPC_DEGREES_PER_RADIAN_R4
822 ierr = gdxyfll(hco%EZscintID, xpos_r4, ypos_r4, &
823 latAdvect_deg_r4, lonAdvect_deg_r4, 1)
824 if (verbose) then
825 write(*,*) 'new xpos_r4 > ni :', &
826 lonIndex0,latIndex0,lonAdvect*MPC_DEGREES_PER_RADIAN_R8, &
827 latAdvect*MPC_DEGREES_PER_RADIAN_R8,stepIndexSF,xpos_r4,ypos_r4
828 end if
829 end if
830
831 ! check if position is west of the grid
832 if (floor(xpos_r4) < 1) then
833 if (verbose) then
834 write(*,*) 'lonIndex0,latIndex0,lon,lat,stepIndexSF,x,y xpos_r4 < 1 :', &
835 lonIndex0,latIndex0,lonAdvect*MPC_DEGREES_PER_RADIAN_R8, &
836 latAdvect*MPC_DEGREES_PER_RADIAN_R8,stepIndexSF,xpos_r4,ypos_r4
837 end if
838 ! subtract 10*epsilon(real*4) to ensure do not go too far due to limited precision
839 lonAdvect = lonAdvect + 2.0D0*MPC_PI_R8 - 10.0D0*real(epsilon(1.0),8)
840 lonAdvect_deg_r4 = real(lonAdvect,4)* MPC_DEGREES_PER_RADIAN_R4
841 latAdvect_deg_r4 = real(latAdvect,4)* MPC_DEGREES_PER_RADIAN_R4
842 ierr = gdxyfll(hco%EZscintID, xpos_r4, ypos_r4, &
843 latAdvect_deg_r4, lonAdvect_deg_r4, 1)
844 if (verbose) then
845 write(*,*) 'new xpos_r4 < 1 :', &
846 lonIndex0,latIndex0,lonAdvect*MPC_DEGREES_PER_RADIAN_R8, &
847 latAdvect*MPC_DEGREES_PER_RADIAN_R8,stepIndexSF,xpos_r4,ypos_r4
848 end if
849 end if
850
851 ! longitude is still outside grid - should not happen!
852 if (floor(xpos_r4) > ni) then
853 write(*,*) '***still outside lonIndex > ni: stepIndexSF,subStepIndex,lonIndex0,latIndex0,x,y,uu=', &
854 stepIndexSF,subStepIndex,lonIndex0,latIndex0,xpos_r4,ypos_r4,uu
855 xpos_r4 = real(ni)
856 lonAdvect = hco%lon(ni)
857 end if
858 if (floor(xpos_r4) < 1) then
859 write(*,*) '***still outside lonIndex < 1 : stepIndexSF,subStepIndex,lonIndex0,latIndex0,x,y,uu=', &
860 stepIndexSF,subStepIndex,lonIndex0,latIndex0,xpos_r4,ypos_r4,uu
861 xpos_r4 = 1.0
862 lonAdvect = hco%lon(1)
863 end if
864
865 ! if position is poleward of last lat circle, ensure valid lat index
866 if (latIndex > nj) then
867 if (verbose) then
868 write(*,*) 'lonIndex0,latIndex0,lon,lat,stepIndexSF,x,y ypos_r4 > nj :', &
869 lonIndex0,latIndex0,lonAdvect*MPC_DEGREES_PER_RADIAN_R8, &
870 latAdvect*MPC_DEGREES_PER_RADIAN_R8,stepIndexSF,xpos_r4,ypos_r4
871 end if
872 ypos_r4 = real(nj)
873 latAdvect = hco%lat(nj)
874 end if
875
876 ! if position is poleward of first lat circle, ensure valid lat index
877 if (latIndex < 1) then
878 if (verbose) then
879 write(*,*) 'lonIndex0,latIndex0,lon,lat,stepIndexSF,x,y ypos_r4 < 1 :', &
880 lonIndex0,latIndex0,lonAdvect*MPC_DEGREES_PER_RADIAN_R8, &
881 latAdvect*MPC_DEGREES_PER_RADIAN_R8,stepIndexSF,xpos_r4,ypos_r4
882 end if
883 ypos_r4 = 1.0
884 latAdvect = hco%lat(1)
885 end if
886
887 ! determine bottom left grid point again after possible adjustments
888 lonIndex = floor(xpos_r4)
889 latIndex = floor(ypos_r4)
890
891 end do ! subStepIndex
892
893 end do ! stepIndexSF
894
895 if (verbose) write(*,*) 'final, initial xpos,ypos', lonIndex0,latIndex0,xpos_r4, ypos_r4
896
897 end SUBROUTINE calcTrajectory
898
899 !--------------------------------------------------------------------------
900 ! calcWeights
901 !--------------------------------------------------------------------------
902 SUBROUTINE calcWeights(lonIndex, latIndex, interpWeight_BL, interpWeight_BR, &
903 interpWeight_TL, interpWeight_TR, xpos_r4, ypos_r4)
904 implicit none
905
906 ! Arguments:
907 integer, intent(out) :: lonIndex
908 integer, intent(out) :: latIndex
909 real(8), intent(out) :: interpWeight_BL
910 real(8), intent(out) :: interpWeight_BR
911 real(8), intent(out) :: interpWeight_TL
912 real(8), intent(out) :: interpWeight_TR
913 real(4), intent(in) :: xpos_r4
914 real(4), intent(in) :: ypos_r4
915
916 ! Locals:
917 real(8) :: delx, dely, sumWeight
918
919 ! Determine bottom left grid point
920 lonIndex = floor(xpos_r4)
921 latIndex = floor(ypos_r4)
922
923 if (lonIndex < 1 .or. lonIndex > hco%ni .or. &
924 latIndex < 1 .or. latIndex > hco%nj ) then
925 write(*,*)
926 write(*,*) 'calcWeights: the input positions are wrong'
927 write(*,*) ' xpos_r4, ypos_r4, lonIndex, latIndex = ', xpos_r4, ypos_r4, lonIndex, latIndex
928 call utl_abort('calcWeights')
929 end if
930
931 ! Compute the surrounding four gridpoint interpolation weights
932 delx = real(xpos_r4,8) - real(lonIndex,8)
933 dely = real(ypos_r4,8) - real(latIndex,8)
934
935 interpWeight_BL = min(max( (1.d0-delx) * (1.d0-dely), 0.0d0), 1.0d0)
936 interpWeight_BR = min(max( delx * (1.d0-dely), 0.0d0), 1.0d0)
937 interpWeight_TL = min(max( (1.d0-delx) * dely , 0.0d0), 1.0d0)
938 interpWeight_TR = min(max( delx * dely , 0.0d0), 1.0d0)
939
940 ! Verification
941 sumWeight = interpWeight_BL + interpWeight_BR + interpWeight_TL + interpWeight_TR
942
943 if (sumWeight > 1.0d0) then
944 write(*,*) 'sumWeight > 1.0 : ', sumWeight
945 write(*,*) ' BL, BR, TL, TR=', &
946 interpWeight_BL, interpWeight_BR, interpWeight_TL, interpWeight_TR
947 write(*,*) ' xpos_r4, ypos_r4, lonIndex, latIndex, delx, dely =', &
948 xpos_r4, ypos_r4, lonIndex, latIndex, delx, dely
949 call utl_abort('calcWeights')
950 end if
951
952 end SUBROUTINE calcWeights
953
954 !--------------------------------------------------------------------------
955 ! adv_ensemble_tl
956 !--------------------------------------------------------------------------
957 SUBROUTINE adv_ensemble_tl( ens, adv, nEns )
958 implicit none
959
960 ! Arguments:
961 type(struct_ens), intent(inout) :: ens
962 type(struct_adv), intent(in) :: adv
963 integer, intent(in) :: nEns
964
965 if ( adv%nLev_M /= ens_getNumLev(ens,'MM') .or. adv%nLev_T /= ens_getNumLev(ens,'TH') ) then
966 call utl_abort('adv_ensemble_tl: vertical levels are not compatible')
967 end if
968
969 if ( ens_getDataKind(ens) == 8 ) then
970 call adv_ensemble_tl_r8( ens, adv, nEns )
971 else if ( ens_getDataKind(ens) == 4 ) then
972 call adv_ensemble_tl_r4( ens, adv, nEns )
973 else
974 call utl_abort('adv_ensemble_tl: ens%dataKind not valid')
975 end if
976
977 END SUBROUTINE adv_ensemble_tl
978
979 !--------------------------------------------------------------------------
980 ! adv_ensemble_tl_r8
981 !--------------------------------------------------------------------------
982 SUBROUTINE adv_ensemble_tl_r8( ens, adv, nEns )
983 implicit none
984
985 ! Arguments:
986 type(struct_ens), intent(inout) :: ens
987 type(struct_adv), intent(in) :: adv
988 integer, intent(in) :: nEns
989
990 ! Locals:
991 real(8), pointer :: ens_oneLev(:,:,:,:)
992 real(8), allocatable :: ens1_mpiglobal_tiles(:,:,:,:)
993 real(8), allocatable :: ens1_mpiglobal(:,:,:)
994 integer :: memberIndex, levIndex, lonIndex, latIndex, kIndex
995 integer :: lonIndex2, latIndex2, lonIndex2_p1, latIndex2_p1, nsize, ierr
996 integer :: procID, procIDx, procIDy, lonIndex_mpiglobal, latIndex_mpiglobal
997 integer :: levTypeIndex, stepIndexAF
998 logical :: gatheringDone
999 character(len=4) :: varName
1000
1001 allocate(ens1_mpiglobal_tiles(nEns,adv%lonPerPE,adv%latPerPE,mmpi_nprocs))
1002 allocate(ens1_mpiglobal(nEns,adv%ni,adv%nj))
1003
1004 do kIndex = 1, ens_getNumK(ens)
1005
1006 levIndex = ens_getLevFromK (ens,kIndex)
1007 varName = ens_getVarNameFromK(ens,kIndex)
1008 if (vnl_varLevelFromVarname(varName) == 'MM') then
1009 levTypeIndex = MMindex
1010 else if (vnl_varLevelFromVarname(varName) == 'TH') then
1011 levTypeIndex = THindex
1012 else if (vnl_varLevelFromVarname(varName) == 'SF') then
1013 levTypeIndex = SFindex
1014 end if
1015
1016 ens_oneLev => ens_getOneLev_r8(ens,kIndex)
1017
1018 gatheringDone = .false.
1019
1020 do stepIndexAF = 1, adv%nTimeStep
1021
1022 if (stepIndexAF == adv%timeStepIndexMainSource) cycle ! no interpolation needed for this time step
1023
1024 if (.not. gatheringDone ) then
1025
1026 ! gather the global field to be interpolated on all tasks
1027 nsize = nEns*adv%lonPerPE*adv%latPerPE
1028 call rpn_comm_allgather(ens_oneLev(1:nEns,adv%timeStepIndexSource(stepIndexAF),:,:), nsize, "mpi_double_precision", &
1029 ens1_mpiglobal_tiles(:,:,:,:), nsize, "mpi_double_precision", &
1030 "GRID", ierr )
1031
1032 ! rearrange gathered fields for convenience
1033 !$OMP PARALLEL DO PRIVATE (procIDy,procIDx,procID,latIndex,lonIndex,latIndex_mpiglobal,lonIndex_mpiglobal,memberIndex)
1034 do procIDy = 0, (mmpi_npey-1)
1035 do procIDx = 0, (mmpi_npex-1)
1036 procID = procIDx + procIDy*mmpi_npex
1037 do latIndex = 1, adv%latPerPE
1038 latIndex_mpiglobal = latIndex + adv%allLatBeg(procIDy+1) - 1
1039 do lonIndex = 1, adv%lonPerPE
1040 lonIndex_mpiglobal = lonIndex + adv%allLonBeg(procIDx+1) - 1
1041 do memberIndex = 1, nEns
1042 ens1_mpiglobal(memberIndex,lonIndex_mpiglobal, latIndex_mpiglobal) = &
1043 ens1_mpiglobal_tiles(memberIndex, lonIndex, latIndex, procID+1)
1044 end do ! memberIndex
1045 end do ! lonIndex
1046 end do ! latIndex
1047 end do ! procIDx
1048 end do ! procIDy
1049 !$OMP END PARALLEL DO
1050
1051 if (adv%singleTimeStepIndexSource) gatheringDone = .true.
1052
1053 end if
1054
1055 !$OMP PARALLEL DO PRIVATE (latIndex,lonIndex,lonIndex2,latIndex2,lonIndex2_p1,latIndex2_p1,memberIndex)
1056 do latIndex = adv%myLatBeg, adv%myLatEnd
1057 do lonIndex = adv%myLonBeg, adv%myLonEnd
1058 lonIndex2 = adv%timeStep(stepIndexAF)%levType(levTypeIndex)%lonIndex(lonIndex,latIndex,levIndex)
1059 latIndex2 = adv%timeStep(stepIndexAF)%levType(levTypeIndex)%latIndex(lonIndex,latIndex,levIndex)
1060 lonIndex2_p1 = mod(lonIndex2,adv%ni)+1 ! assume periodic
1061 latIndex2_p1 = min(latIndex2+1,adv%nj)
1062 do memberIndex = 1, nEns
1063 ens_oneLev(memberIndex,stepIndexAF,lonIndex,latIndex) = &
1064 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BL(lonIndex,latIndex,levIndex)* &
1065 ens1_mpiglobal(memberIndex, lonIndex2 ,latIndex2 ) + &
1066 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BR(lonIndex,latIndex,levIndex)* &
1067 ens1_mpiglobal(memberIndex, lonIndex2_p1,latIndex2 ) + &
1068 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TL(lonIndex,latIndex,levIndex)* &
1069 ens1_mpiglobal(memberIndex, lonIndex2 ,latIndex2_p1) + &
1070 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TR(lonIndex,latIndex,levIndex)* &
1071 ens1_mpiglobal(memberIndex, lonIndex2_p1,latIndex2_p1)
1072 end do ! memberIndex
1073 end do ! lonIndex
1074 end do ! latIndex
1075 !$OMP END PARALLEL DO
1076
1077 end do ! stepIndexAF
1078
1079 end do ! kIndex
1080
1081 deallocate(ens1_mpiglobal_tiles)
1082 deallocate(ens1_mpiglobal)
1083
1084 END SUBROUTINE adv_ensemble_tl_r8
1085
1086 !--------------------------------------------------------------------------
1087 ! adv_ensemble_tl_r8
1088 !--------------------------------------------------------------------------
1089 SUBROUTINE adv_ensemble_tl_r4( ens, adv, nEns )
1090 implicit none
1091
1092 ! Arguments:
1093 type(struct_ens), intent(inout) :: ens
1094 type(struct_adv), intent(in) :: adv
1095 integer, intent(in) :: nEns
1096
1097 ! Locals:
1098 real(4), pointer :: ens_oneLev(:,:,:,:)
1099 real(4), allocatable :: ens1_mpiglobal_tiles(:,:,:,:)
1100 real(4), allocatable :: ens1_mpiglobal(:,:,:)
1101 integer :: memberIndex, levIndex, lonIndex, latIndex, kIndex
1102 integer :: lonIndex2, latIndex2, lonIndex2_p1, latIndex2_p1, nsize, ierr
1103 integer :: procID, procIDx, procIDy, lonIndex_mpiglobal, latIndex_mpiglobal
1104 integer :: levTypeIndex, stepIndexAF
1105 logical :: gatheringDone
1106 character(len=4) :: varName
1107
1108 allocate(ens1_mpiglobal_tiles(nEns,adv%lonPerPE,adv%latPerPE,mmpi_nprocs))
1109 allocate(ens1_mpiglobal(nEns,adv%ni,adv%nj))
1110
1111 do kIndex = 1, ens_getNumK(ens)
1112
1113 levIndex = ens_getLevFromK (ens,kIndex)
1114 varName = ens_getVarNameFromK(ens,kIndex)
1115 if (vnl_varLevelFromVarname(varName) == 'MM') then
1116 levTypeIndex = MMindex
1117 else if (vnl_varLevelFromVarname(varName) == 'TH') then
1118 levTypeIndex = THindex
1119 else if (vnl_varLevelFromVarname(varName) == 'SF') then
1120 levTypeIndex = SFindex
1121 end if
1122
1123 ens_oneLev => ens_getOneLev_r4(ens,kIndex)
1124
1125 gatheringDone = .false.
1126
1127 do stepIndexAF = 1, adv%nTimeStep
1128
1129 if (stepIndexAF == adv%timeStepIndexMainSource) cycle ! no interpolation needed for this time step
1130
1131 if (.not. gatheringDone ) then
1132
1133 ! gather the global field to be interpolated on all tasks
1134 nsize = nEns*adv%lonPerPE*adv%latPerPE
1135 call rpn_comm_allgather(ens_oneLev(1:nEns,adv%timeStepIndexSource(stepIndexAF),:,:), nsize, "mpi_real4", &
1136 ens1_mpiglobal_tiles(:,:,:,:), nsize, "mpi_real4", &
1137 "GRID", ierr )
1138
1139 ! rearrange gathered fields for convenience
1140 !$OMP PARALLEL DO PRIVATE (procIDy,procIDx,procID,latIndex,lonIndex,latIndex_mpiglobal,lonIndex_mpiglobal,memberIndex)
1141 do procIDy = 0, (mmpi_npey-1)
1142 do procIDx = 0, (mmpi_npex-1)
1143 procID = procIDx + procIDy*mmpi_npex
1144 do latIndex = 1, adv%latPerPE
1145 latIndex_mpiglobal = latIndex + adv%allLatBeg(procIDy+1) - 1
1146 do lonIndex = 1, adv%lonPerPE
1147 lonIndex_mpiglobal = lonIndex + adv%allLonBeg(procIDx+1) - 1
1148 do memberIndex = 1, nEns
1149 ens1_mpiglobal(memberIndex,lonIndex_mpiglobal, latIndex_mpiglobal) = &
1150 ens1_mpiglobal_tiles(memberIndex, lonIndex, latIndex, procID+1)
1151 end do ! memberIndex
1152 end do ! lonIndex
1153 end do ! latIndex
1154 end do ! procIDx
1155 end do ! procIDy
1156 !$OMP END PARALLEL DO
1157
1158 if (adv%singleTimeStepIndexSource) gatheringDone = .true.
1159
1160 end if
1161
1162 !$OMP PARALLEL DO PRIVATE (latIndex,lonIndex,lonIndex2,latIndex2,lonIndex2_p1,latIndex2_p1,memberIndex)
1163 do latIndex = adv%myLatBeg, adv%myLatEnd
1164 do lonIndex = adv%myLonBeg, adv%myLonEnd
1165 lonIndex2 = adv%timeStep(stepIndexAF)%levType(levTypeIndex)%lonIndex(lonIndex,latIndex,levIndex)
1166 latIndex2 = adv%timeStep(stepIndexAF)%levType(levTypeIndex)%latIndex(lonIndex,latIndex,levIndex)
1167 lonIndex2_p1 = mod(lonIndex2,adv%ni)+1 ! assume periodic
1168 latIndex2_p1 = min(latIndex2+1,adv%nj)
1169 do memberIndex = 1, nEns
1170 ens_oneLev(memberIndex,stepIndexAF,lonIndex,latIndex) = &
1171 real(adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BL(lonIndex,latIndex,levIndex),4)* &
1172 ens1_mpiglobal(memberIndex, lonIndex2 ,latIndex2 ) + &
1173 real(adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BR(lonIndex,latIndex,levIndex),4)* &
1174 ens1_mpiglobal(memberIndex, lonIndex2_p1,latIndex2 ) + &
1175 real(adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TL(lonIndex,latIndex,levIndex),4)* &
1176 ens1_mpiglobal(memberIndex, lonIndex2 ,latIndex2_p1) + &
1177 real(adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TR(lonIndex,latIndex,levIndex),4)* &
1178 ens1_mpiglobal(memberIndex, lonIndex2_p1,latIndex2_p1)
1179 end do ! memberIndex
1180 end do ! lonIndex
1181 end do ! latIndex
1182 !$OMP END PARALLEL DO
1183
1184 end do ! stepIndexAF
1185
1186 end do ! kIndex
1187
1188 deallocate(ens1_mpiglobal_tiles)
1189 deallocate(ens1_mpiglobal)
1190
1191 END SUBROUTINE adv_ensemble_tl_r4
1192
1193 !--------------------------------------------------------------------------
1194 ! adv_ensemble_ad
1195 !--------------------------------------------------------------------------
1196 SUBROUTINE adv_ensemble_ad( ens, adv, nEns )
1197 implicit none
1198
1199 ! Arguments:
1200 type(struct_ens), intent(inout) :: ens
1201 type(struct_adv), intent(in) :: adv
1202 integer, intent(in) :: nEns
1203
1204 ! Locals:
1205 real(8), pointer :: ens_oneLev(:,:,:,:)
1206 real(8), allocatable :: ens1_mpiglobal(:,:,:)
1207 real(8), allocatable :: ens1_mpiglobal_tiles(:,:,:,:)
1208 real(8), allocatable :: ens1_mpiglobal_tiles2(:,:,:,:)
1209 integer :: memberIndex, levIndex, lonIndex, latIndex, kIndex
1210 integer :: lonIndex2, latIndex2, lonIndex2_p1, latIndex2_p1, nsize, ierr
1211 integer :: procID, procIDx, procIDy, lonIndex_mpiglobal, latIndex_mpiglobal
1212 integer :: levTypeIndex, stepIndexAF
1213 character(len=4) :: varName
1214
1215 if ( .not. adv%singleTimeStepIndexSource ) then
1216 call utl_abort('adv_ensemble_ad cannot deal with multiple timeStep index source')
1217 end if
1218 if ( ens_getDataKind(ens) /= 8 ) then
1219 call utl_abort('adv_ensemble_ad can only deal with double precision (real8) ensembleStateVector')
1220 end if
1221 if ( adv%nLev_M /= ens_getNumLev(ens,'MM') .or. adv%nLev_T /= ens_getNumLev(ens,'TH') ) then
1222 call utl_abort('adv_ensemble_ad: vertical levels are not compatible')
1223 end if
1224 if ( .not. nlat_equalAcrossMpiTasks .or. .not. nlon_equalAcrossMpiTasks) then
1225 call utl_abort('adv_ensemble_ad can only deal with even nlon and nlat across all MPI tasks')
1226 end if
1227
1228 allocate(ens1_mpiglobal(nEns,adv%ni,adv%nj))
1229 allocate(ens1_mpiglobal_tiles (nEns,adv%lonPerPE,adv%latPerPE,mmpi_nprocs))
1230 allocate(ens1_mpiglobal_tiles2(nEns,adv%lonPerPE,adv%latPerPE,mmpi_nprocs))
1231
1232 do kIndex = 1, ens_getNumK(ens)
1233
1234 levIndex = ens_getLevFromK (ens,kIndex)
1235 varName = ens_getVarNameFromK(ens,kIndex)
1236 if (vnl_varLevelFromVarname(varName) == 'MM') then
1237 levTypeIndex = MMindex
1238 else if (vnl_varLevelFromVarname(varName) == 'TH') then
1239 levTypeIndex = THindex
1240 else if (vnl_varLevelFromVarname(varName) == 'SF') then
1241 levTypeIndex = SFindex
1242 end if
1243
1244 ens1_mpiglobal(:,:,:) = 0.0d0
1245 ens_oneLev => ens_getOneLev_r8(ens,kIndex)
1246
1247 do latIndex = adv%myLatBeg, adv%myLatEnd
1248 do lonIndex = adv%myLonBeg, adv%myLonEnd
1249 do stepIndexAF = 1, adv%nTimeStep
1250 if (stepIndexAF == adv%timeStepIndexMainSource) cycle ! no interpolation needed for this time step
1251 ! this is the bottom-left grid point
1252 lonIndex2 = adv%timeStep(stepIndexAF)%levType(levTypeIndex)%lonIndex(lonIndex,latIndex,levIndex)
1253 latIndex2 = adv%timeStep(stepIndexAF)%levType(levTypeIndex)%latIndex(lonIndex,latIndex,levIndex)
1254 lonIndex2_p1 = mod(lonIndex2,adv%ni)+1 ! assume periodic
1255 latIndex2_p1 = min(latIndex2+1,adv%nj)
1256 !$OMP PARALLEL DO PRIVATE(memberIndex)
1257 do memberIndex = 1, nEns
1258 ens1_mpiglobal(memberIndex, lonIndex2 ,latIndex2) = &
1259 ens1_mpiglobal(memberIndex, lonIndex2 ,latIndex2) + &
1260 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BL(lonIndex,latIndex,levIndex)* &
1261 ens_oneLev(memberIndex,stepIndexAF,lonIndex,latIndex)
1262 ens1_mpiglobal(memberIndex, lonIndex2_p1,latIndex2) = &
1263 ens1_mpiglobal(memberIndex, lonIndex2_p1,latIndex2) + &
1264 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BR(lonIndex,latIndex,levIndex)* &
1265 ens_oneLev(memberIndex,stepIndexAF,lonIndex,latIndex)
1266 ens1_mpiglobal(memberIndex, lonIndex2 ,latIndex2_p1) = &
1267 ens1_mpiglobal(memberIndex, lonIndex2 ,latIndex2_p1) + &
1268 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TL(lonIndex,latIndex,levIndex)* &
1269 ens_oneLev(memberIndex,stepIndexAF,lonIndex,latIndex)
1270 ens1_mpiglobal(memberIndex, lonIndex2_p1,latIndex2_p1) = &
1271 ens1_mpiglobal(memberIndex, lonIndex2_p1,latIndex2_p1) + &
1272 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TR(lonIndex,latIndex,levIndex)* &
1273 ens_oneLev(memberIndex,stepIndexAF,lonIndex,latIndex)
1274 end do ! memberIndex
1275 !$OMP END PARALLEL DO
1276 end do ! stepIndexAF
1277 end do ! lonIndex
1278 end do ! latIndex
1279
1280 ! redistribute the global initial time field across mpi tasks by tiles
1281 !$OMP PARALLEL DO PRIVATE(procIDy,procIDx,procID,latIndex,latIndex_mpiglobal,lonIndex,lonIndex_mpiglobal,memberIndex)
1282 do procIDy = 0, (mmpi_npey-1)
1283 do procIDx = 0, (mmpi_npex-1)
1284 procID = procIDx + procIDy*mmpi_npex
1285 do latIndex = 1, adv%latPerPE
1286 latIndex_mpiglobal = latIndex + adv%allLatBeg(procIDy+1) - 1
1287 do lonIndex = 1, adv%lonPerPE
1288 lonIndex_mpiglobal = lonIndex + adv%allLonBeg(procIDx+1) - 1
1289 do memberIndex = 1, nEns
1290 ens1_mpiglobal_tiles(memberIndex, lonIndex, latIndex, procID+1) = &
1291 ens1_mpiglobal(memberIndex, lonIndex_mpiglobal, latIndex_mpiglobal)
1292 end do ! memberIndex
1293 end do ! lonIndex
1294 end do ! latIndex
1295 end do ! procIDx
1296 end do ! procIDy
1297 !$OMP END PARALLEL DO
1298
1299 nsize = nEns*adv%lonPerPE*adv%latPerPE
1300 if (mmpi_nprocs > 1) then
1301 call rpn_comm_alltoall(ens1_mpiglobal_tiles, nsize,"mpi_double_precision", &
1302 ens1_mpiglobal_tiles2,nsize,"mpi_double_precision","GRID",ierr)
1303 else
1304 ens1_mpiglobal_tiles2(:,:,:,1) = ens1_mpiglobal_tiles(:,:,:,1)
1305 end if
1306
1307 do procID = 0, (mmpi_nprocs-1)
1308 !$OMP PARALLEL DO PRIVATE(latIndex,latIndex2,lonIndex,lonIndex2,memberIndex)
1309 do latIndex = 1, adv%latPerPE
1310 latIndex2= latIndex + adv%myLatBeg - 1
1311 do lonIndex = 1, adv%lonPerPE
1312 lonIndex2 = lonIndex + adv%myLonBeg - 1
1313 do memberIndex = 1, nEns
1314 ens_oneLev(memberIndex, adv%timeStepIndexMainSource, lonIndex2, latIndex2) = &
1315 ens_oneLev(memberIndex, adv%timeStepIndexMainSource, lonIndex2, latIndex2) + &
1316 ens1_mpiglobal_tiles2(memberIndex, lonIndex, latIndex, procID+1)
1317 end do ! memberIndex
1318 end do ! lonIndex
1319 end do ! latIndex
1320 !$OMP END PARALLEL DO
1321 end do ! procID
1322
1323 end do ! levIndex
1324
1325 deallocate(ens1_mpiglobal)
1326 deallocate(ens1_mpiglobal_tiles)
1327 deallocate(ens1_mpiglobal_tiles2)
1328
1329 END SUBROUTINE adv_ensemble_ad
1330
1331 !--------------------------------------------------------------------------
1332 ! adv_statevector_tl
1333 !--------------------------------------------------------------------------
1334 SUBROUTINE adv_statevector_tl( statevector, adv)
1335 implicit none
1336
1337 ! Arguments:
1338 type(struct_gsv), intent(inout) :: statevector
1339 type(struct_adv), intent(in) :: adv
1340
1341 ! Locals:
1342 real(8), pointer :: field4D(:,:,:,:)
1343 real(8), allocatable :: field2D_mpiglobal_tiles(:,:,:)
1344 real(8), allocatable :: field2D_mpiglobal(:,:)
1345 integer :: levIndex, lonIndex, latIndex, kIndex
1346 integer :: lonIndex2, latIndex2, lonIndex2_p1, latIndex2_p1, nsize, ierr
1347 integer :: procID, procIDx, procIDy, lonIndex_mpiglobal, latIndex_mpiglobal
1348 integer :: levTypeIndex, stepIndexAF
1349 logical :: gatheringDone
1350 character(len=4) :: varName
1351
1352 if ( gsv_getDataKind(statevector) /= 8 ) then
1353 call utl_abort('adv_statevector_tl can only deal with double precision (real8) gridStateVector')
1354 end if
1355 if ( adv%nLev_M /= statevector%vco%nLev_M .or. adv%nLev_T /= statevector%vco%nLev_T ) then
1356 call utl_abort('adv_statevector_tl: vertical levels are not compatible')
1357 end if
1358
1359 call utl_tmg_start(100,'--ADV_GSV')
1360
1361 allocate(field2D_mpiglobal_tiles(adv%lonPerPE,adv%latPerPE,mmpi_nprocs))
1362 allocate(field2D_mpiglobal(adv%ni,adv%nj))
1363
1364 do kIndex = 1, gsv_getNumK(statevector)
1365
1366 levIndex = gsv_getLevFromK (statevector,kIndex)
1367 varName = gsv_getVarNameFromK(statevector,kIndex)
1368 if (vnl_varLevelFromVarname(varName) == 'MM') then
1369 levTypeIndex = MMindex
1370 else if (vnl_varLevelFromVarname(varName) == 'TH') then
1371 levTypeIndex = THindex
1372 else if (vnl_varLevelFromVarname(varName) == 'SF') then
1373 levTypeIndex = SFindex
1374 end if
1375
1376 call gsv_getField(statevector, field4D, varName)
1377
1378 gatheringDone = .false.
1379
1380 do stepIndexAF = 1, adv%nTimeStep
1381
1382 if (stepIndexAF == adv%timeStepIndexMainSource) cycle ! no interpolation needed for this time step
1383
1384 if (.not. gatheringDone ) then
1385
1386 ! gather the global field to be interpolated on all tasks
1387 call rpn_comm_barrier('GRID',ierr)
1388 call utl_tmg_start(101,'----ADV_GSV_Comm')
1389 nsize = adv%lonPerPE*adv%latPerPE
1390 call rpn_comm_allgather(field4D(:,:,levIndex,adv%timeStepIndexSource(stepIndexAF)), nsize, "mpi_double_precision", &
1391 field2D_mpiglobal_tiles(:,:,:), nsize, "mpi_double_precision", &
1392 "GRID", ierr )
1393 call utl_tmg_stop(101)
1394
1395 ! rearrange gathered fields for convenience
1396 call utl_tmg_start(102,'----ADV_GSV_Shuffling')
1397 !$OMP PARALLEL DO PRIVATE (procIDy,procIDx,procID,latIndex,lonIndex,latIndex_mpiglobal,lonIndex_mpiglobal)
1398 do procIDy = 0, (mmpi_npey-1)
1399 do procIDx = 0, (mmpi_npex-1)
1400 procID = procIDx + procIDy*mmpi_npex
1401 do latIndex = 1, adv%latPerPE
1402 latIndex_mpiglobal = latIndex + adv%allLatBeg(procIDy+1) - 1
1403 do lonIndex = 1, adv%lonPerPE
1404 lonIndex_mpiglobal = lonIndex + adv%allLonBeg(procIDx+1) - 1
1405 field2D_mpiglobal(lonIndex_mpiglobal, latIndex_mpiglobal) = &
1406 field2D_mpiglobal_tiles(lonIndex, latIndex, procID+1)
1407 end do ! lonIndex
1408 end do ! latIndex
1409 end do ! procIDx
1410 end do ! procIDy
1411 !$OMP END PARALLEL DO
1412 call utl_tmg_stop(102)
1413
1414 if (adv%singleTimeStepIndexSource) gatheringDone = .true.
1415
1416 end if
1417
1418 call utl_tmg_start(103,'----ADV_GSV_Calc')
1419
1420 !$OMP PARALLEL DO PRIVATE (latIndex,lonIndex,lonIndex2,latIndex2,lonIndex2_p1,latIndex2_p1)
1421 do latIndex = adv%myLatBeg, adv%myLatEnd
1422 do lonIndex = adv%myLonBeg, adv%myLonEnd
1423 lonIndex2 = adv%timeStep(stepIndexAF)%levType(levTypeIndex)%lonIndex(lonIndex,latIndex,levIndex)
1424 latIndex2 = adv%timeStep(stepIndexAF)%levType(levTypeIndex)%latIndex(lonIndex,latIndex,levIndex)
1425 lonIndex2_p1 = mod(lonIndex2,adv%ni)+1 ! assume periodic
1426 latIndex2_p1 = min(latIndex2+1,adv%nj)
1427 field4D(lonIndex,latIndex,levIndex,stepIndexAF) = &
1428 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BL(lonIndex,latIndex,levIndex)* &
1429 field2D_mpiglobal(lonIndex2 ,latIndex2 ) + &
1430 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BR(lonIndex,latIndex,levIndex)* &
1431 field2D_mpiglobal(lonIndex2_p1,latIndex2 ) + &
1432 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TL(lonIndex,latIndex,levIndex)* &
1433 field2D_mpiglobal(lonIndex2 ,latIndex2_p1) + &
1434 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TR(lonIndex,latIndex,levIndex)* &
1435 field2D_mpiglobal(lonIndex2_p1,latIndex2_p1)
1436 end do ! lonIndex
1437 end do ! latIndex
1438 !$OMP END PARALLEL DO
1439
1440 call utl_tmg_stop(103)
1441
1442 end do ! stepIndexAF
1443
1444 end do ! kIndex
1445
1446 deallocate(field2D_mpiglobal_tiles)
1447 deallocate(field2D_mpiglobal)
1448
1449 call utl_tmg_stop(100)
1450
1451 END SUBROUTINE adv_statevector_tl
1452
1453 !--------------------------------------------------------------------------
1454 ! adv_statevector_ad
1455 !--------------------------------------------------------------------------
1456 SUBROUTINE adv_statevector_ad( statevector, adv)
1457 implicit none
1458
1459 ! Arguments:
1460 type(struct_gsv), intent(inout) :: statevector
1461 type(struct_adv), intent(in) :: adv
1462
1463 ! Locals:
1464 real(8), pointer :: field4D(:,:,:,:)
1465 real(8), allocatable :: field2D_mpiglobal(:,:)
1466 real(8), allocatable :: field2D_mpiglobal_tiles (:,:,:)
1467 real(8), allocatable :: field2D_mpiglobal_tiles2(:,:,:)
1468 integer :: levIndex, lonIndex, latIndex, kIndex
1469 integer :: lonIndex2, latIndex2, lonIndex2_p1, latIndex2_p1, nsize, ierr
1470 integer :: procID, procIDx, procIDy, lonIndex_mpiglobal, latIndex_mpiglobal
1471 integer :: levTypeIndex, stepIndexAF
1472 character(len=4) :: varName
1473
1474 if ( adv%singleTimeStepIndexSource ) then
1475 call utl_abort('adv_statevector_ad cannot work for singleTimeStepIndexSource')
1476 end if
1477 if ( gsv_getDataKind(statevector) /= 8 ) then
1478 call utl_abort('adv_statevector_ad can only deal with double precision (real8) ensembleStateVector')
1479 end if
1480 if ( adv%nLev_M /= statevector%vco%nLev_M .or. adv%nLev_T /= statevector%vco%nLev_T ) then
1481 call utl_abort('adv_statevector_ad: vertical levels are not compatible')
1482 end if
1483 if ( .not. nlat_equalAcrossMpiTasks .or. .not. nlon_equalAcrossMpiTasks) then
1484 call utl_abort('adv_ensemble_ad can only deal with even nlon and nlat across all MPI tasks')
1485 end if
1486
1487 allocate(field2D_mpiglobal(adv%ni,adv%nj))
1488 allocate(field2D_mpiglobal_tiles (adv%lonPerPE,adv%latPerPE,mmpi_nprocs))
1489 allocate(field2D_mpiglobal_tiles2(adv%lonPerPE,adv%latPerPE,mmpi_nprocs))
1490
1491 do kIndex = 1, gsv_getNumK(statevector)
1492
1493 levIndex = gsv_getLevFromK (statevector,kIndex)
1494 varName = gsv_getVarNameFromK(statevector,kIndex)
1495 if (vnl_varLevelFromVarname(varName) == 'MM') then
1496 levTypeIndex = MMindex
1497 else if (vnl_varLevelFromVarname(varName) == 'TH') then
1498 levTypeIndex = THindex
1499 else if (vnl_varLevelFromVarname(varName) == 'SF') then
1500 levTypeIndex = SFindex
1501 end if
1502
1503 call gsv_getField(statevector, field4D, varName)
1504
1505 do stepIndexAF = 1, adv%nTimeStep
1506
1507 if (stepIndexAF == adv%timeStepIndexMainSource) cycle ! no interpolation needed for this time step
1508
1509 field2D_mpiglobal(:,:) = 0.0d0
1510
1511 do latIndex = adv%myLatBeg, adv%myLatEnd
1512 do lonIndex = adv%myLonBeg, adv%myLonEnd
1513 ! this is the bottom-left grid point
1514 lonIndex2 = adv%timeStep(stepIndexAF)%levType(levTypeIndex)%lonIndex(lonIndex,latIndex,levIndex)
1515 latIndex2 = adv%timeStep(stepIndexAF)%levType(levTypeIndex)%latIndex(lonIndex,latIndex,levIndex)
1516 lonIndex2_p1 = mod(lonIndex2,adv%ni)+1 ! assume periodic
1517 latIndex2_p1 = min(latIndex2+1,adv%nj)
1518 field2D_mpiglobal(lonIndex2 ,latIndex2) = &
1519 field2D_mpiglobal(lonIndex2 ,latIndex2) + &
1520 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BL(lonIndex,latIndex,levIndex)* &
1521 field4D(lonIndex,latIndex,levIndex,stepIndexAF)
1522 field2D_mpiglobal(lonIndex2_p1,latIndex2) = &
1523 field2D_mpiglobal(lonIndex2_p1,latIndex2) + &
1524 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_BR(lonIndex,latIndex,levIndex)* &
1525 field4D(lonIndex,latIndex,levIndex,stepIndexAF)
1526 field2D_mpiglobal(lonIndex2 ,latIndex2_p1) = &
1527 field2D_mpiglobal(lonIndex2 ,latIndex2_p1) + &
1528 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TL(lonIndex,latIndex,levIndex)* &
1529 field4D(lonIndex,latIndex,levIndex,stepIndexAF)
1530 field2D_mpiglobal(lonIndex2_p1,latIndex2_p1) = &
1531 field2D_mpiglobal(lonIndex2_p1,latIndex2_p1) + &
1532 adv%timeStep(stepIndexAF)%levType(levTypeIndex)%interpWeight_TR(lonIndex,latIndex,levIndex)* &
1533 field4D(lonIndex,latIndex,levIndex,stepIndexAF)
1534 end do ! lonIndex
1535 end do ! latIndex
1536
1537 ! redistribute the global initial time field across mpi tasks by tiles
1538 !$OMP PARALLEL DO PRIVATE(procIDy,procIDx,procID,latIndex,latIndex_mpiglobal,lonIndex,lonIndex_mpiglobal)
1539 do procIDy = 0, (mmpi_npey-1)
1540 do procIDx = 0, (mmpi_npex-1)
1541 procID = procIDx + procIDy*mmpi_npex
1542 do latIndex = 1, adv%latPerPE
1543 latIndex_mpiglobal = latIndex + adv%allLatBeg(procIDy+1) - 1
1544 do lonIndex = 1, adv%lonPerPE
1545 lonIndex_mpiglobal = lonIndex + adv%allLonBeg(procIDx+1) - 1
1546 field2D_mpiglobal_tiles(lonIndex, latIndex, procID+1) = &
1547 field2D_mpiglobal(lonIndex_mpiglobal, latIndex_mpiglobal)
1548 end do ! lonIndex
1549 end do ! latIndex
1550 end do ! procIDx
1551 end do ! procIDy
1552 !$OMP END PARALLEL DO
1553
1554 nsize = adv%lonPerPE*adv%latPerPE
1555 if (mmpi_nprocs > 1) then
1556 call rpn_comm_alltoall(field2D_mpiglobal_tiles, nsize,"mpi_double_precision", &
1557 field2D_mpiglobal_tiles2,nsize,"mpi_double_precision","GRID",ierr)
1558 else
1559 field2D_mpiglobal_tiles2(:,:,1) = field2D_mpiglobal_tiles(:,:,1)
1560 end if
1561
1562 field4D(:, :, levIndex, adv%timeStepIndexSource(stepIndexAF)) = 0.d0
1563
1564 do procID = 0, (mmpi_nprocs-1)
1565 !$OMP PARALLEL DO PRIVATE(latIndex,latIndex2,lonIndex,lonIndex2)
1566 do latIndex = 1, adv%latPerPE
1567 latIndex2= latIndex + adv%myLatBeg - 1
1568 do lonIndex = 1, adv%lonPerPE
1569 lonIndex2 = lonIndex + adv%myLonBeg - 1
1570 field4D(lonIndex2, latIndex2, levIndex, adv%timeStepIndexSource(stepIndexAF)) = &
1571 field4D(lonIndex2, latIndex2, levIndex, adv%timeStepIndexSource(stepIndexAF)) + &
1572 field2D_mpiglobal_tiles2(lonIndex, latIndex, procID+1)
1573 end do ! lonIndex
1574 end do ! latIndex
1575 !$OMP END PARALLEL DO
1576 end do ! procID
1577
1578 end do ! stepIndexAF
1579
1580 end do ! levIndex
1581
1582 deallocate(field2D_mpiglobal)
1583 deallocate(field2D_mpiglobal_tiles)
1584 deallocate(field2D_mpiglobal_tiles2)
1585
1586 END SUBROUTINE adv_statevector_ad
1587
1588END MODULE advection_mod
