1module stateToColumn_mod
2 ! MODULE stateToColumn_mod (prefix='s2c' category='4. Data Object transformations')
3 !
4 !:Purpose: Non-linear, tangent-linear and adjoint versions of
5 ! horizontal-temporal interpolation between a gridStateVector object
6 ! and a columnData object.
7 !
8 use mathPhysConstants_mod
9 use earthConstants_mod
10 use mpi, only : mpi_status_size ! this is the mpi library module
11 use midasMpi_mod
12 use codePrecision_mod
13 use gridStateVector_mod
14 use obsSpaceData_mod
15 use columnData_mod
16 use horizontalCoord_mod
17 use obsTimeInterp_mod
18 use windRotation_mod
19 use utilities_mod
20 use gridVariableTransforms_mod
21 use varNameList_mod
22 use slantProfileLatLon_mod
23 use tovsNL_mod
24 use codtyp_mod
25 use getGridPosition_mod
26 use kdTree2_mod
27 use calcHeightAndPressure_mod
28 use humidityLimits_mod
29
30 implicit none
31 save
32 private
33
34 ! public routines
35 public :: s2c_tl, s2c_ad, s2c_nl
36 public :: s2c_bgcheck_bilin, s2c_getFootprintRadius, s2c_getWeightsAndGridPointIndexes
37 public :: s2c_deallocInterpInfo
38
39 ! private module variables and derived types
40
41 type struct_stepProcData
42 ! lat-lon location of observations to be interpolated
43 real(8), pointer :: allLat(:,:) => null() ! (headerUsed, kIndex)
44 real(8), pointer :: allLon(:,:) => null() ! (headerUsed, kIndex)
45 ! lat-lon location on rotated grid of observations to be interpolated
46 real(8), pointer :: allLatRot(:,:,:) => null() ! (subGrid, headerUsed, kIndex)
47 real(8), pointer :: allLonRot(:,:,:) => null() ! (subGrid, headerUsed, kIndex)
48 ! actual headerIndex, since the headerUsed is only for those obs with a non-zero interp weight
49 integer, pointer :: allHeaderIndex(:) => null() ! (headerUsed)
50 ! depotIndexBeg/End contain first/last indices into depots of interpolation weights and lat/lon indices
51 integer, pointer :: depotIndexBeg(:,:,:) => null() ! (subGrid, headerUsed, kIndex)
52 integer, pointer :: depotIndexEnd(:,:,:) => null() ! (subGrid, headerUsed, kIndex)
53 end type struct_stepProcData
54
55 type struct_interpInfo
56 logical :: initialized = .false.
57 type(struct_hco), pointer :: hco => null() ! horizontal grid object
58 type(struct_uvr), pointer :: uvr => null() ! windRotation object
59 type(struct_oti), pointer :: oti => null() ! obsTimeInterp object
60
61 ! number of obs headers on each proc having a non-zero interp weight for each stepIndex (headerUsed)
62 integer, pointer :: allNumHeaderUsed(:,:) => null() ! (step, proc)
63
64 ! structure containing all interpolation information that depends on (proc, step)
65 type(struct_stepProcData), allocatable :: stepProcData(:,:) ! (proc, step)
66
67 ! interpolation weights and lat/lon indices are accessed via the 'stepProcData%depotIndexBeg/End'
68 real(8), allocatable :: interpWeightDepot(:) ! (depotIndex)
69 integer, pointer :: latIndexDepot(:) ! (depotIndex)
70 integer, pointer :: lonIndexDepot(:) ! (depotIndex)
71 character(len=2) :: inputStateVectorType
72 end type struct_interpInfo
73
74 type(struct_interpInfo), target :: interpInfo_tlad, interpInfo_nl
75 type(kdtree2), pointer :: tree_nl => null()
76 type(kdtree2), pointer :: tree_tlad => null()
77
78 character(len=20), parameter :: timeInterpType_tlad = 'LINEAR' ! hardcoded type of time interpolation for increment
79
80 integer, parameter :: maxNumWrites = 50
81 logical, parameter :: verbose = .false.
82
83 ! "special" values of the footprint radius
84 real(4), parameter :: nearestNeighbourFootprint = -3.0
85 real(4), parameter :: lakeFootprint = -2.0
86 real(4), parameter :: bilinearFootprint = -1.0
87 integer, parameter :: maxNumLocalGridptsSearch = 1000
88 integer, parameter :: minNumLocalGridptsSearch = 8
89
90 ! namelist variables
91 logical :: slantPath_TO_nl ! choose to use slant path for non-linear radiance operator
92 logical :: slantPath_TO_tlad ! choose to use slant path for linearized radiance operators
93 logical :: slantPath_RO_nl ! choose to use slant path for non-linear GPS-RO operator
94 logical :: slantPath_RA_nl ! choose to use slant path for non-linear radar operator
95 logical :: calcHeightPressIncrOnColumn ! choose to compute Z and P increment on column, instead of grid
96 logical :: useFootprintForTovs ! choose to use a horizontal footprint operator for radiance obs
97 logical :: rejectObsNonMonotonicPressure ! choose to reject obs when interpolated column pressure is non-monotonic
98
99 integer, external :: get_max_rss
100
101contains
102
103
104 !---------------------------------------------------------
105 ! pressureProfileMonotonicityCheck
106 !---------------------------------------------------------
107 subroutine pressureProfileMonotonicityCheck(obsSpaceData, column)
108 !
109 !:Purpose: Check for non monotonic pressure profiles that can be computed in slantpathmode
110 !
111 implicit none
112
113 ! Arguments:
114 type(struct_obs), intent(inout) :: obsSpaceData
115 type(struct_columnData), intent(inout) :: column
116
117 ! Locals:
118 integer, parameter :: numWriteMax = 10
119 integer :: headerIndex, bodyIndex, iterationCount, singularIndex, levelIndex
120 integer :: pressureVarIndex
121 integer :: nlv
122 integer :: numWrites
123 real(8), pointer :: pressureProfile(:)
124 logical :: monotonicProfile
125 integer, parameter :: nPressureVar =2
126 character(len=4), parameter :: pressureVarList(nPressureVar)=['P_T ', 'P_M ']
127
128 write(*,*) ' '
129 write(*,*) 'pressureProfileMonotonicityCheck: START'
130 write(*,*) ' '
131 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
132
133 numWrites = 0
134
135 call obs_set_current_header_list(obsSpaceData, 'TO')
136 HEADER: do
137 headerIndex = obs_getHeaderIndex(obsSpaceData)
138 if (headerIndex < 0) exit HEADER
139 do pressureVarIndex = 1, nPressureVar
140 pressureProfile => col_getColumn(column, headerIndex, pressureVarList(pressureVarIndex))
141 nlv = size(pressureProfile)
142 monotonicProfile = .true.
143 iterationCount = 0
144 iterationLoop:do
145 singularIndex = -1
146 levelSearch:do levelIndex = 1, nlv - 1
147 if ( pressureProfile(levelIndex) > pressureProfile(levelIndex+1)) then
148 singularIndex = levelIndex
149 exit levelSearch
150 end if
151 end do levelSearch
152 if ( singularIndex == -1 ) exit iterationLoop !regular profile or correction OK
153 iterationCount = iterationCount + 1
154 if (iterationCount == 1) then
155 monotonicProfile = .false.
156 if (numWrites < numWriteMax) then
157 numWrites = numWrites + 1
158 write(*,*) 'pressureProfileMonotonicityCheck: found non monotonic pressure profile:', &
159 pressureVarList(pressureVarIndex), pressureProfile
160 end if
161 end if
162 if (singularIndex == 1) then !should never happen
163 write(*,*) 'pressureProfileMonotonicityCheck: ', pressureProfile(1:2)
164 call utl_abort('pressureProfileMonotonicityCheck: profile in the wrong order ?' &
165 // pressureVarList(pressureVarIndex))
166 end if
167 pressureProfile(singularIndex) = 0.5d0 * ( pressureProfile(singularIndex - 1) + pressureProfile(singularIndex + 1) )
168 write(*,*) 'pressureProfileMonotonicityCheck: profile iteration', &
169 pressureVarList(pressureVarIndex), iterationCount
170 end do iterationLoop
171
172 ! if requested reject the corrected profile
173 if (.not. monotonicProfile .and. rejectObsNonMonotonicPressure) then
174 call obs_headSet_i(obsSpaceData,OBS_ST1,headerIndex, &
175 ibset( obs_headElem_i(obsSpaceData,OBS_ST1,headerIndex), 05))
176 call obs_set_current_body_list(obsSpaceData, headerIndex)
177 BODY: do
178 bodyIndex = obs_getBodyIndex(obsSpaceData)
179 if (bodyIndex < 0) exit BODY
180 if (rejectObsNonMonotonicPressure) then
181 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
182 call obs_bodySet_i(obsSpaceData, OBS_FLG, bodyIndex, &
183 ibset(obs_bodyElem_i(obsSpaceData, OBS_FLG, bodyIndex),9))
184 end if
185 end do BODY
186 end if
187 end do ! loop on pressure variables
188
189 end do HEADER
190
191 write(*,*) 'pressureProfileMonotonicityCheck: END'
192
193 end subroutine pressureProfileMonotonicityCheck
194
195 !---------------------------------------------------------
196 ! latlonChecksAnlGrid
197 !---------------------------------------------------------
198 subroutine latlonChecksAnlGrid(obsSpaceData, hco_core, moveObsAtPole)
199 !
200 !:Purpose: Check the lat/lon of observations and modify if necessary
201 !
202 implicit none
203
204 ! Arguments:
205 type(struct_obs) , intent(inout) :: obsSpaceData
206 type(struct_hco), pointer, intent(in) :: hco_core
207 logical , intent(in) :: moveObsAtPole
208
209 ! Locals:
210 integer :: headerIndex, ierr
211 integer :: idata, idatend, jdata, subGridIndex
212 real(4) :: lat_r4, lon_r4, lat_deg_r4, lon_deg_r4
213 real(4) :: xpos_r4, ypos_r4, xpos2_r4, ypos2_r4
214 real(4) :: xposLowerBoundAnl_r4, xposUpperBoundAnl_r4
215 real(8) :: lat_r8, lon_r8
216 integer, save :: numWrites = 0
217 integer :: gdllfxy
218
219 write(*,*) ' '
220 write(*,*) 'latlonChecksAnlGrid: STARTING'
221 write(*,*) ' '
222 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
223
224 !
225 !- Get the Analysis Grid structure
226 !
227 if ( hco_core % global ) then
228 xposLowerBoundAnl_r4 = - huge(1.0) ! no limit since grid is global (periodic)
229 xposUpperBoundAnl_r4 = + huge(1.0) ! no limit since grid is global (periodic)
230 else
231 xposLowerBoundAnl_r4 = 1.0
232 xposUpperBoundAnl_r4 = real(hco_core % ni)
233 end if
234
235 header_loop: do headerIndex=1, obs_numheader(obsSpaceData)
236
237 !- Get LatLon of observation location
238 lat_r8 = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
239 lon_r8 = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
240 lat_r4 = real(lat_r8,4)
241 lon_r4 = real(lon_r8,4)
242 if (lon_r4.lt.0.0 ) lon_r4 = lon_r4 + 2.0*MPC_PI_R4
243 if (lon_r4.ge.2.*MPC_PI_R4) lon_r4 = lon_r4 - 2.0*MPC_PI_R4
244
245 lat_deg_r4 = lat_r4 * MPC_DEGREES_PER_RADIAN_R4 ! Radian To Degree
246 lon_deg_r4 = lon_r4 * MPC_DEGREES_PER_RADIAN_R4
247
248 !
249 !- Find the position in the analysis grid
250 !
251 ierr = gpos_getPositionXY( hco_core % EZscintID, &
252 xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
253 lat_deg_r4, lon_deg_r4, subGridIndex )
254
255 !- Test if the obs is outside the analysis grid
256 if ( xpos_r4 < xposLowerBoundAnl_r4 .or. &
257 xpos_r4 > xposUpperBoundAnl_r4 .or. &
258 ypos_r4 < 1.0 .or. &
259 ypos_r4 > real(hco_core % nj) ) then
260
261 if ( hco_core % global ) then
262
263 if ( moveObsAtPole ) then
264 ! Modify latitude if we have an observation at or near the poles
265 write(*,*) ''
266 write(*,*) 'latlonChecksAnlGrid: Moving OBS inside the GLOBAL ANALYSIS grid, ', headerIndex
267 write(*,*) ' true position : ', lat_deg_r4, lon_deg_r4, ypos_r4, xpos_r4
268
269 !- Move the observation to the nearest grid point
270 if ( ypos_r4 < 1.0 ) ypos_r4 = 1.0
271 if ( ypos_r4 > real(hco_core % nj) ) ypos_r4 = real(hco_core % nj)
272
273 ierr = gdllfxy( hco_core % EZscintID, & ! IN
274 lat_deg_r4, lon_deg_r4, & ! OUT
275 xpos_r4, ypos_r4, 1) ! IN
276
277 write(*,*) ' new position : ', lat_deg_r4, lon_deg_r4, ypos_r4, xpos_r4
278
279 lat_r8 = real(lat_deg_r4,8) * MPC_RADIANS_PER_DEGREE_R8
280 lon_r8 = real(lon_deg_r4,8) * MPC_RADIANS_PER_DEGREE_R8
281 call obs_headSet_r(obsSpaceData,OBS_LAT,headerIndex, lat_r8) ! IN
282 call obs_headSet_r(obsSpaceData,OBS_LON,headerIndex, lon_r8) ! IN
283 else
284 write(*,*)
285 write(*,*) 'latlonChecksAnlGrid: OBS outside the GLOBAL ANALYSIS grid, but NOT moved, ', headerIndex
286 end if
287
288 else
289 ! The observation is outside the domain
290 ! In LAM Analysis mode we must discard this observation
291 numWrites = numWrites + 1
292 if (numWrites < maxNumWrites) then
293 write(*,*) 'latlonChecksAnlGrid: Rejecting OBS outside the LAM ANALYSIS grid domain, ', headerIndex
294 write(*,*) ' position : ', lat_deg_r4, lon_deg_r4, ypos_r4, xpos_r4
295 else if (numWrites == maxNumWrites) then
296 write(*,*) 'latlonChecksAnlGrid: More rejects, but reached maximum number of writes to the listing.'
297 end if
298
299 idata = obs_headElem_i(obsSpaceData,OBS_RLN,headerIndex)
300 idatend = obs_headElem_i(obsSpaceData,OBS_NLV,headerIndex) + idata -1
301 do jdata = idata, idatend
302 call obs_bodySet_i(obsSpaceData,OBS_ASS,JDATA, obs_notAssimilated)
303 end do
304 call obs_headSet_i(obsSpaceData,OBS_ST1,headerIndex, &
305 ibset( obs_headElem_i(obsSpaceData,OBS_ST1,headerIndex), 05))
306 end if
307
308 end if
309
310 end do header_loop
311
312 write(*,*) 'latlonChecksAnlGrid: END'
313
314 end subroutine latlonChecksAnlGrid
315
316 !---------------------------------------------------------
317 ! s2c_setupInterpInfo
318 !---------------------------------------------------------
319 subroutine s2c_setupInterpInfo( interpInfo, obsSpaceData, stateVector, &
320 headerIndexBeg, headerIndexEnd, &
321 timeInterpType, rejectOutsideObs, &
322 inputStateVectorType, lastCall_opt )
323 !:Purpose: Setup all of the information needed to quickly
324 ! perform the horizontal interpolation to the observation
325 ! locations.
326 !
327 implicit none
328
329 ! Arguments:
330 type(struct_interpInfo), intent(out) :: interpInfo
331 type(struct_obs) , intent(inout) :: obsSpaceData
332 type(struct_gsv), target, intent(in) :: stateVector
333 integer , intent(in) :: headerIndexBeg
334 integer , intent(in) :: headerIndexEnd
335 logical , intent(in) :: rejectOutsideObs
336 character(len=*) , intent(in) :: timeInterpType
337 character(len=*) , intent(in) :: inputStateVectorType
338 logical, optional , intent(in) :: lastCall_opt
339
340 ! Locals:
341 type(struct_gsv) :: stateVector_VarsLevs_1Step, stateVector_Tiles_allVar_1Step
342 type(struct_gsv) :: stateVector_Tiles_1Step
343 type(struct_gsv), save :: stateVector_1Step
344 type(struct_gsv), pointer :: stateVector_Tiles_ptr
345 integer :: numHeader, numHeaderUsedMax, headerIndex, headerUsedIndex
346 integer :: kIndex, kIndexCount, myKBeg
347 integer :: numStep, stepIndex, fnom, fclos, nulnam, ierr
348 integer :: procIndex, niP1, numGridptTotal, numHeaderUsed
349 integer :: subGridIndex, subGridForInterp, numSubGridsForInterp
350 real(8) :: latRot, lonRot, lat, lon
351 real(4) :: lon_r4, lat_r4, lon_deg_r4, lat_deg_r4
352 real(4) :: xpos_r4, ypos_r4, xpos2_r4, ypos2_r4
353 real(4) :: footprintRadius_r4 ! (metres)
354 integer, allocatable :: numGridpt(:), allNumHeaderUsed(:,:)
355 integer, allocatable :: allHeaderIndex(:,:,:), headerIndexVec(:,:)
356 real(8), allocatable :: lat_send_r8(:,:), lat_recv_r8(:,:), lon_send_r8(:,:), lon_recv_r8(:,:)
357 real(4), allocatable :: footprintRadiusVec_r4(:), allFootprintRadius_r4(:,:,:)
358 real(8), allocatable :: allLatOneLev(:,:)
359 real(8), allocatable :: allLonOneLev(:,:)
360 character(len=4), pointer :: varNames(:)
361 character(len=4) :: varLevel
362 real(8), allocatable :: latColumn(:,:), lonColumn(:,:)
363 real(8), allocatable :: latLev_T(:), lonLev_T(:), latLev_M(:), lonLev_M(:)
364 real(8) :: latLev_S, lonLev_S
365 real(4), pointer :: height3D_r4_ptr1(:,:,:), height3D_r4_ptr2(:,:,:)
366 real(4), save, pointer :: height3D_T_r4(:,:,:), height3D_M_r4(:,:,:)
367 real(8), pointer :: height3D_r8_ptr1(:,:,:)
368 real(kdkind), allocatable :: positionArray(:,:)
369 integer :: sendsizes(mmpi_nprocs), recvsizes(mmpi_nprocs), senddispls(mmpi_nprocs)
370 integer :: recvdispls(mmpi_nprocs), allkBeg(mmpi_nprocs)
371 integer :: codeType, nlev_T, nlev_M, levIndex
372 integer :: lonIndex, latIndex, gridIndex
373 integer :: maxkcount, numkToSend, numTovsUsingFootprint, numAllTovs
374 logical :: doSlantPath, SlantTO, SlantRO, SlantRA, firstHeaderSlantPathTO, firstHeaderSlantPathRO, firstHeaderSlantPathRA
375 logical :: doSetup3dHeights, lastCall
376 logical, save :: nmlAlreadyRead = .false.
377 type(kdtree2), pointer :: tree
378
379 namelist /nams2c/ slantPath_TO_nl, slantPath_TO_tlad, slantPath_RO_nl, slantPath_RA_nl, calcHeightPressIncrOnColumn
380 namelist /nams2c/ useFootprintForTovs, rejectObsNonMonotonicPressure
381
382 write(*,*) 's2c_setupInterpInfo: STARTING'
383 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
384
385 write(*,*) 's2c_setupInterpInfo: inputStateVectorType=', inputStateVectorType
386
387 if ( present(lastCall_opt) ) then
388 lastCall = lastCall_opt
389 else
390 lastCall = .false.
391 end if
392
393 if ( .not. nmlAlreadyRead ) then
394 nmlAlreadyRead = .true.
395
396 ! default values
397 slantPath_TO_nl = .false.
398 slantPath_TO_tlad = .false.
399 slantPath_RO_nl = .false.
400 slantPath_RA_nl = .false.
401 calcHeightPressIncrOnColumn = .false.
402 useFootprintForTovs = .false.
403 rejectObsNonMonotonicPressure =.true.
404
405 if ( .not. utl_isNamelistPresent('NAMS2C','./flnml') ) then
406 if ( mmpi_myid == 0 ) then
407 write(*,*) 's2c_setupInterpInfo: nams2c is missing in the namelist.'
408 write(*,*) ' The default values will be taken.'
409 end if
410
411 else
412 ! reading namelist variables
413 nulnam = 0
414 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
415 read(nulnam, nml = nams2c, iostat = ierr)
416 if ( ierr /= 0 ) call utl_abort('s2c_setupInterpInfo: Error reading namelist')
417 ierr = fclos(nulnam)
418 end if
419 if ( mmpi_myid == 0 ) write(*, nml = nams2c)
420 end if
421
422 doSlantPath = .false.
423 SlantTO = .false.
424 SlantRO = .false.
425 if ( slantPath_TO_nl .and. inputStateVectorType == 'nl' ) then
426 doSlantPath = .true.
427 SlantTO = .true.
428 endif
429 if ( slantPath_TO_tlad .and. inputStateVectorType /= 'nl' ) then
430 doSlantPath = .true.
431 SlantTO = .true.
432 endif
433 if ( slantPath_RO_nl .and. inputStateVectorType == 'nl' ) then
434 doSlantPath = .true.
435 SlantRO = .true.
436 endif
437 if ( slantPath_RA_nl .and. inputStateVectorType == 'nl' ) then
438 doSlantPath = .true.
439 SlantRA = .true.
440 endif
441 write(*,*) 's2c_setupInterpInfo: doSlantPath, SlantTO, SlantRO, SlantRA = ', &
442 doSlantPath, SlantTO, SlantRO, SlantRA
443
444 numStep = stateVector%numStep
445 numHeader = headerIndexEnd - headerIndexBeg + 1
446
447 call oti_setup(interpInfo%oti, obsSpaceData, numStep, &
448 headerIndexBeg, headerIndexEnd, &
449 interpType_opt=timeInterpType, flagObsOutside_opt=.true.)
450
451 if ((stateVector%heightSfcPresent) .and. ( mmpi_myid == 0)) then
452 mykBeg = 0
453 else
454 mykBeg = stateVector%mykBeg
455 end if
456
457 call rpn_comm_allgather(mykBeg, 1,'mpi_integer', &
458 allkBeg,1,'mpi_integer','grid',ierr)
459
460 ! Allow for periodicity in Longitude for global Gaussian grid
461 if ( stateVector%hco%grtyp == 'G' .or. &
462 (stateVector%hco%grtyp == 'Z' .and. stateVector%hco%global) ) then
463 niP1 = stateVector%ni + 1
464 else
465 niP1 = stateVector%ni
466 end if
467
468 ! First count the number of headers for each stepIndex
469 allocate(allNumHeaderUsed(numStep,mmpi_nprocs))
470 allNumHeaderUsed(:,:) = 0
471 do stepIndex = 1, numStep
472 numHeaderUsed = 0
473
474 header_loop1: do headerIndex = headerIndexBeg, headerIndexEnd
475
476 ! if obs inside window, but zero weight for current stepIndex then skip it
477 if ( oti_getTimeInterpWeight(interpInfo%oti,headerIndex,stepIndex) == 0.0d0 ) then
478 cycle header_loop1
479 end if
480
481 numHeaderUsed = numHeaderUsed + 1
482
483 end do header_loop1
484 ! gather the number of obs over all processors for each timestep
485 call rpn_comm_allgather(numHeaderUsed, 1, 'MPI_INTEGER', &
486 allNumHeaderUsed(stepIndex,:), 1, 'MPI_INTEGER', &
487 'GRID',ierr)
488
489 end do
490
491 numHeaderUsedMax = maxval(allNumHeaderUsed(:,:))
492 write(*,*) 's2c_setupInterpInfo: numHeaderUsedMax = ', numHeaderUsedMax
493
494 ! temporary arrays
495 allocate(headerIndexVec(numHeaderUsedMax,numStep))
496 allocate(footprintRadiusVec_r4(numHeaderUsedMax))
497 headerIndexVec(:,:) = 0
498
499 ! copy the horizontal grid object
500 interpInfo%hco => stateVector%hco
501
502 ! setup the information for wind rotation
503 if ( (gsv_varExist(varName='UU') .or. gsv_varExist(varName='VV')) .and. &
504 stateVector%hco%rotated ) then
505 call uvr_Setup( interpInfo%uvr, & ! INOUT
506 stateVector%hco ) ! IN
507 end if
508
509 allocate(interpInfo%stepProcData(mmpi_nprocs,numStep))
510 do stepIndex = 1,numStep
511 do procIndex = 1, mmpi_nprocs
512 allocate(interpInfo%stepProcData(procIndex,stepIndex)%allLat(allNumHeaderUsed(stepIndex,procIndex),mykBeg:stateVector%mykEnd))
513 allocate(interpInfo%stepProcData(procIndex,stepIndex)%allLon(allNumHeaderUsed(stepIndex,procIndex),mykBeg:stateVector%mykEnd))
514 interpInfo%stepProcData(procIndex,stepIndex)%allLat(:,:) = 0.0d0
515 interpInfo%stepProcData(procIndex,stepIndex)%allLon(:,:) = 0.0d0
516
517 allocate(interpInfo%stepProcData(procIndex,stepIndex)%allHeaderIndex(allNumHeaderUsed(stepIndex,procIndex)))
518 interpInfo%stepProcData(procIndex,stepIndex)%allHeaderIndex(:) = 0
519
520 allocate(interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(interpInfo%hco%numSubGrid,numHeaderUsedMax,mykBeg:stateVector%mykEnd))
521 allocate(interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(interpInfo%hco%numSubGrid,numHeaderUsedMax,mykBeg:stateVector%mykEnd))
522 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(:,:,:) = 0
523 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(:,:,:) = -1
524 end do
525 end do
526
527 ! allocate arrays that will be returned
528 allocate(interpInfo%allNumHeaderUsed(numStep,mmpi_nprocs))
529 allocate(allLatOneLev(numHeaderUsedMax,mmpi_nprocs))
530 allocate(allLonOneLev(numHeaderUsedMax,mmpi_nprocs))
531 allocate(allFootprintRadius_r4(numHeaderUsedMax,numStep,mmpi_nprocs))
532 allocate(numGridpt(interpInfo%hco%numSubGrid))
533 allFootprintRadius_r4(:,:,:) = bilinearFootprint
534 interpInfo%allNumHeaderUsed(:,:) = allNumHeaderUsed(:,:)
535
536 if ( interpInfo%hco%rotated ) then
537 do stepIndex = 1, numStep
538 do procIndex = 1, mmpi_nprocs
539 allocate(interpInfo%stepProcData(procIndex,stepIndex)%allLatRot(interpInfo%hco%numSubGrid,allNumHeaderUsed(stepIndex,procIndex),mykBeg:stateVector%mykEnd))
540 allocate(interpInfo%stepProcData(procIndex,stepIndex)%allLonRot(interpInfo%hco%numSubGrid,allNumHeaderUsed(stepIndex,procIndex),mykBeg:stateVector%mykEnd))
541 interpInfo%stepProcData(procIndex,stepIndex)%allLatRot(:,:,:) = 0.0d0
542 interpInfo%stepProcData(procIndex,stepIndex)%allLonRot(:,:,:) = 0.0d0
543 end do
544 end do
545 end if
546
547 nlev_T = gsv_getNumLev(stateVector,'TH')
548 nlev_M = gsv_getNumLev(stateVector,'MM')
549
550 doSetup3dHeights = doSlantPath .and. &
551 .not. gsv_isAllocated(stateVector_1Step) .and. &
552 gsv_varExist(stateVector,'Z_T') .and. &
553 gsv_varExist(stateVector,'Z_M')
554
555 ! prepare for extracting the 3D height for slant-path calculation
556 if ( doSetup3dHeights ) then
557
558 write(*,*) 's2c_setupInterpInfo: extracting 3D heights for slant-path for ', inputStateVectorType
559
560 if ( inputStateVectorType == 'nl' ) then
561 nullify(varNames)
562 call gsv_varNamesList(varNames, stateVector)
563 call gsv_allocate( stateVector_VarsLevs_1Step, 1, &
564 stateVector%hco, stateVector%vco, &
565 mpi_local_opt=.true., mpi_distribution_opt='VarsLevs', &
566 dataKind_opt=4, varNames_opt=varNames )
567
568 call gsv_getField(stateVector,height3D_r4_ptr1)
569 call gsv_getField(stateVector_VarsLevs_1Step,height3D_r4_ptr2)
570 height3D_r4_ptr2(:,:,:) = height3D_r4_ptr1(:,:,:)
571
572 call gsv_allocate( stateVector_Tiles_allVar_1Step, 1, &
573 stateVector%hco, stateVector%vco, &
574 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
575 dataKind_opt=4, varNames_opt=varNames )
576
577 call gsv_transposeVarsLevsToTiles( stateVector_VarsLevs_1Step, stateVector_Tiles_allVar_1Step )
578 call gsv_deallocate(statevector_VarsLevs_1Step)
579
580 call gsv_allocate( stateVector_Tiles_1Step, 1, &
581 stateVector%hco, stateVector%vco, &
582 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
583 dataKind_opt=4, varNames_opt=(/'Z_M','Z_T'/) )
584
585 call gsv_getField(stateVector_Tiles_allVar_1Step,height3D_r4_ptr1,'Z_T')
586 call gsv_getField(stateVector_Tiles_1Step,height3D_r4_ptr2,'Z_T')
587 height3D_r4_ptr2(:,:,:) = height3D_r4_ptr1(:,:,:)
588
589 call gsv_getField(stateVector_Tiles_allVar_1Step,height3D_r4_ptr1,'Z_M')
590 call gsv_getField(stateVector_Tiles_1Step,height3D_r4_ptr2,'Z_M')
591 height3D_r4_ptr2(:,:,:) = height3D_r4_ptr1(:,:,:)
592
593 call gsv_deallocate(stateVector_Tiles_allVar_1Step)
594
595 else
596 stateVector_Tiles_ptr => gvt_getStateVectorTrial('height')
597
598 call gsv_allocate( stateVector_Tiles_1Step, 1, &
599 stateVector%hco, stateVector%vco, &
600 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
601 dataKind_opt=4, varNames_opt=(/'Z_M','Z_T'/) )
602
603 call gsv_getField(stateVector_Tiles_ptr,height3D_r8_ptr1,'Z_T')
604 call gsv_getField(stateVector_Tiles_1Step,height3D_r4_ptr2,'Z_T')
605 height3D_r4_ptr2(:,:,:) = height3D_r8_ptr1(:,:,:)
606
607 call gsv_getField(stateVector_Tiles_ptr,height3D_r8_ptr1,'Z_M')
608 call gsv_getField(stateVector_Tiles_1Step,height3D_r4_ptr2,'Z_M')
609 height3D_r4_ptr2(:,:,:) = height3D_r8_ptr1(:,:,:)
610
611 end if ! inputStateVectorType
612
613 ! Communicate 3D height fields onto all mpi tasks
614 call gsv_allocate( stateVector_1Step, 1, &
615 stateVector%hco, stateVector%vco, &
616 mpi_local_opt=.false., &
617 dataKind_opt=4, varNames_opt=(/'Z_M','Z_T'/) )
618 call utl_tmg_start(32,'------s2c_Slant')
619 call gsv_transposeTilesToMpiGlobal(stateVector_1Step, stateVector_Tiles_1Step)
620 call utl_tmg_stop(32)
621 call gsv_getField(stateVector_1Step,height3D_T_r4,'Z_T')
622 call gsv_getField(stateVector_1Step,height3D_M_r4,'Z_M')
623
624 write(*,*) 's2c_setupInterpInfo, height3D_T_r4='
625 write(*,*) height3D_T_r4(1,1,:)
626 write(*,*) 's2c_setupInterpInfo, height3D_M_r4='
627 write(*,*) height3D_M_r4(1,1,:)
628
629 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
630 end if ! doSlantPath
631
632 ! get observation lat-lon and footprint radius onto all mpi tasks
633 step_loop2: do stepIndex = 1, numStep
634 numHeaderUsed = 0
635
636 footprintRadiusVec_r4(:) = bilinearFootprint
637
638 header_loop2: do headerIndex = headerIndexBeg, headerIndexEnd
639
640 ! if obs inside window, but zero weight for current stepIndex then skip it
641 if ( oti_getTimeInterpWeight(interpInfo%oti, headerIndex, stepIndex) == 0.0d0 ) then
642 cycle header_loop2
643 end if
644
645 numHeaderUsed = numHeaderUsed + 1
646 headerIndexVec(numHeaderUsed,stepIndex) = headerIndex
647
648 footprintRadiusVec_r4(numHeaderUsed) = s2c_getFootprintRadius(obsSpaceData, &
649 stateVector, headerIndex)
650
651 end do header_loop2
652
653 call rpn_comm_allgather(footprintRadiusVec_r4, numHeaderUsedMax, 'MPI_REAL4', &
654 allFootprintRadius_r4(:,stepIndex,:), numHeaderUsedMax, 'MPI_REAL4', &
655 'GRID', ierr)
656
657 allocate(latColumn(numHeaderUsedMax,allkBeg(1):stateVector%nk))
658 allocate(lonColumn(numHeaderUsedMax,allkBeg(1):stateVector%nk))
659 latColumn(:,:) = 0.0d0
660 lonColumn(:,:) = 0.0d0
661
662 if ( doSlantPath .and. &
663 gsv_varExist(stateVector,'Z_T') .and. &
664 gsv_varExist(stateVector,'Z_M') ) then
665
666 allocate(latLev_T(nlev_T))
667 allocate(lonLev_T(nlev_T))
668 allocate(latLev_M(nlev_M))
669 allocate(lonLev_M(nlev_M))
670 latLev_T(:) = 0.0d0
671 lonLev_T(:) = 0.0d0
672 latLev_M(:) = 0.0d0
673 lonLev_M(:) = 0.0d0
674
675 firstHeaderSlantPathTO = .true.
676 firstHeaderSlantPathRO = .true.
677 firstHeaderSlantPathRA = .true.
678 header_loop3: do headerUsedIndex = 1, numHeaderUsed
679 headerIndex = headerIndexVec(headerUsedIndex,stepIndex)
680
681 !- Get LatLon of observation location
682 lat_r4 = real(obs_headElem_r(obsSpaceData, OBS_LAT, headerIndex), 4)
683 lon_r4 = real(obs_headElem_r(obsSpaceData, OBS_LON, headerIndex), 4)
684 if (lon_r4 < 0.0 ) lon_r4 = lon_r4 + 2.0*MPC_PI_R4
685 if (lon_r4 >= 2.0*MPC_PI_R4) lon_r4 = lon_r4 - 2.0*MPC_PI_R4
686
687 codeType = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
688
689 if ( tvs_isIdBurpTovs(codeType) .and. SlantTO ) then
690 if ( firstHeaderSlantPathTO ) then
691 write(*,'(a,i3,a,i8)') 's2c_setupInterpInfo: start slant-path for TOVS. stepIndex = ', &
692 stepIndex,' and numHeaderUsed = ',numHeaderUsed
693 firstHeaderSlantPathTO = .false.
694 end if
695
696 ! calculate lat/lon along the line of sight
697 call utl_tmg_start(32,'------s2c_Slant')
698 call slp_calcLatLonTovs( obsSpaceData, stateVector%hco, headerIndex, & ! IN
699 height3D_T_r4, height3D_M_r4, & ! IN
700 latLev_T, lonLev_T, & ! OUT
701 latLev_M, lonLev_M, & ! OUT
702 latLev_S, lonLev_S ) ! OUT
703 call utl_tmg_stop(32)
704
705 else if (codeType == codtyp_get_codtyp('ro') .and. SlantRO ) then
706 if ( firstHeaderSlantPathRO ) then
707 write(*,'(a,i3,a,i8)') 's2c_setupInterpInfo: start slant-path for RO. stepIndex = ', &
708 stepIndex,' and numHeaderUsed = ',numHeaderUsed
709 firstHeaderSlantPathRO = .false.
710 end if
711
712 ! Calculate lat/lon along the GPSRO obs
713 call utl_tmg_start(32,'------s2c_Slant')
714 call slp_calcLatLonRO( obsSpaceData, stateVector%hco, headerIndex, & ! IN
715 height3D_T_r4, height3D_M_r4, & ! IN
716 latLev_T, lonLev_T, & ! OUT
717 latLev_M, lonLev_M, & ! OUT
718 latLev_S, lonLev_S ) ! OUT
719 call utl_tmg_stop(32)
720 else if (codeType == codtyp_get_codtyp('radar') .and. SlantRA ) then
721 if ( firstHeaderSlantPathRA ) then
722 write(*,'(a,i3,a,i8)') 's2c_setupInterpInfo: start slant-path for RADAR. stepIndex=', &
723 stepIndex,' and numHeaderUsed=',numHeaderUsed
724 firstHeaderSlantPathRA = .false.
725 end if
726
727 ! calculate lat/lon along the radar beam obs
728 call slp_calcLatLonRadar( obsSpaceData, stateVector%hco, headerIndex, & ! IN
729 height3D_T_r4, height3D_M_r4, & ! IN
730 latLev_T, lonLev_T, & ! OUT
731 latLev_M, lonLev_M, & ! OUT
732 latLev_S, lonLev_S ) ! OUT
733 else
734
735 latLev_T(:) = real(lat_r4,8)
736 lonLev_T(:) = real(lon_r4,8)
737 latLev_M(:) = real(lat_r4,8)
738 lonLev_M(:) = real(lon_r4,8)
739 latLev_S = real(lat_r4,8)
740 lonLev_S = real(lon_r4,8)
741 end if
742
743 ! check if the slanted lat/lon is inside the domain
744 call latlonChecks ( obsSpaceData, stateVector%hco, & ! IN
745 headerIndex, rejectOutsideObs, & ! IN
746 latLev_T, lonLev_T, & ! IN/OUT
747 latLev_M, lonLev_M, & ! IN/OUT
748 latLev_S, lonLev_S ) ! IN/OUT
749
750 ! put the lat/lon from TH/MM levels to kIndex
751 do kIndex = allkBeg(1), stateVector%nk
752 if ( kIndex == 0 ) then
753 varLevel = 'SF'
754 else
755 levIndex = gsv_getLevFromK(stateVector,kIndex)
756 varLevel = vnl_varLevelFromVarname(gsv_getVarNameFromK(stateVector,kIndex))
757 end if
758
759 if ( varLevel == 'TH' ) then
760 latColumn(headerUsedIndex,kIndex) = latLev_T(levIndex)
761 lonColumn(headerUsedIndex,kIndex) = lonLev_T(levIndex)
762 else if ( varLevel == 'MM' ) then
763 latColumn(headerUsedIndex,kIndex) = latLev_M(levIndex)
764 lonColumn(headerUsedIndex,kIndex) = lonLev_M(levIndex)
765 else if ( varLevel == 'SF' ) then
766 latColumn(headerUsedIndex,kIndex) = latLev_S
767 lonColumn(headerUsedIndex,kIndex) = lonLev_S
768 else
769 call utl_abort('s2c_setupInterpInfo: unknown value of varLevel')
770 end if
771
772 end do
773
774 end do header_loop3
775
776 ! MPI communication for the slant-path lat/lon
777
778 maxkCount = maxval(stateVector%allkCount(:) + stateVector%allkBeg(:) - allkBeg(:))
779 numkToSend = min(mmpi_nprocs,stateVector%nk)
780
781 allocate(lat_recv_r8(numHeaderUsedMax,mmpi_nprocs))
782 lat_recv_r8(:,:) = 0.0d0
783 allocate(lat_send_r8(numHeaderUsedMax,mmpi_nprocs))
784 lat_send_r8(:,:) = 0.0d0
785 allocate(lon_recv_r8(numHeaderUsedMax,mmpi_nprocs))
786 lon_recv_r8(:,:) = 0.0d0
787 allocate(lon_send_r8(numHeaderUsedMax,mmpi_nprocs))
788 lon_send_r8(:,:) = 0.0d0
789
790 ! only send the data from tasks with data, same amount to all
791 sendsizes(:) = 0
792 do procIndex = 1, numkToSend
793 sendsizes(procIndex) = numHeaderUsed
794 end do
795 senddispls(1) = 0
796 do procIndex = 2, mmpi_nprocs
797 senddispls(procIndex) = senddispls(procIndex-1) + numHeaderUsedMax
798 end do
799
800 recvdispls(1) = 0
801 do procIndex = 2, mmpi_nprocs
802 recvdispls(procIndex) = recvdispls(procIndex-1) + numHeaderUsedMax
803 end do
804
805 ! loop to send (at most) 1 level to (at most) all other mpi tasks
806 do kIndexCount = 1, maxkCount
807
808 sendsizes(:) = 0
809 do procIndex = 1, mmpi_nprocs
810 ! compute kIndex value being sent
811 kIndex = kIndexCount + allkBeg(procIndex) - 1
812 if ( kIndex <= stateVector%allkEnd(procIndex) ) then
813 if( procIndex > numkToSend ) then
814 write(*,*) 'procIndex, numkToSend = ', procIndex, numkToSend
815 call utl_abort('ERROR: with numkToSend?')
816 end if
817
818 lat_send_r8(1:numHeaderUsed,procIndex) = latColumn(1:numHeaderUsed,kIndex)
819 lon_send_r8(1:numHeaderUsed,procIndex) = lonColumn(1:numHeaderUsed,kIndex)
820 sendsizes(procIndex) = numHeaderUsed
821 else
822 sendsizes(procIndex) = 0
823 end if
824 end do
825
826 ! all tasks recv only from those with data
827 kIndex = kIndexCount + mykBeg - 1
828 if ( kIndex <= stateVector%mykEnd ) then
829 do procIndex = 1, mmpi_nprocs
830 recvsizes(procIndex) = allNumHeaderUsed(stepIndex,procIndex)
831 end do
832 else
833 recvsizes(:) = 0
834 end if
835
836 call mpi_alltoallv(lat_send_r8, sendsizes, senddispls, mmpi_datyp_real8, &
837 lat_recv_r8, recvsizes, recvdispls, mmpi_datyp_real8, &
838 mmpi_comm_grid, ierr)
839 call mpi_alltoallv(lon_send_r8, sendsizes, senddispls, mmpi_datyp_real8, &
840 lon_recv_r8, recvsizes, recvdispls, mmpi_datyp_real8, &
841 mmpi_comm_grid, ierr)
842
843 do procIndex = 1, mmpi_nprocs
844 ! all tasks copy the received step data into correct slot
845 kIndex = kIndexCount + mykBeg - 1
846 if ( kIndex <= stateVector%mykEnd ) then
847 interpInfo%stepProcData(procIndex,stepIndex)%allLat(:,kIndex) = &
848 lat_recv_r8(1:allNumHeaderUsed(stepIndex,procIndex),procIndex)
849 interpInfo%stepProcData(procIndex,stepIndex)%allLon(:,kIndex) = &
850 lon_recv_r8(1:allNumHeaderUsed(stepIndex,procIndex),procIndex)
851 end if
852 end do
853
854 end do ! kIndexCount
855
856 deallocate(lon_send_r8)
857 deallocate(lon_recv_r8)
858 deallocate(lat_send_r8)
859 deallocate(lat_recv_r8)
860
861 deallocate(latLev_T)
862 deallocate(lonLev_T)
863 deallocate(latLev_M)
864 deallocate(lonLev_M)
865
866 else ! not doSlantPath
867
868 allocate(latLev_T(1))
869 allocate(lonLev_T(1))
870 allocate(latLev_M(1))
871 allocate(lonLev_M(1))
872 latLev_T(:) = 0.0d0
873 lonLev_T(:) = 0.0d0
874 latLev_M(:) = 0.0d0
875 lonLev_M(:) = 0.0d0
876
877 do headerUsedIndex = 1, numHeaderUsed
878 headerIndex = headerIndexVec(headerUsedIndex,stepIndex)
879
880 !- Get LatLon of observation location
881 lat_r4 = real(obs_headElem_r(obsSpaceData, OBS_LAT, headerIndex), 4)
882 lon_r4 = real(obs_headElem_r(obsSpaceData, OBS_LON, headerIndex), 4)
883 if (lon_r4 < 0.0 ) lon_r4 = lon_r4 + 2.0*MPC_PI_R4
884 if (lon_r4 >= 2.0*MPC_PI_R4) lon_r4 = lon_r4 - 2.0*MPC_PI_R4
885
886 latLev_T(:) = real(lat_r4,8)
887 lonLev_T(:) = real(lon_r4,8)
888 latLev_M(:) = real(lat_r4,8)
889 lonLev_M(:) = real(lon_r4,8)
890
891 ! check if the lat/lon is inside the domain
892 call latlonChecks ( obsSpaceData, stateVector%hco, & ! IN
893 headerIndex, rejectOutsideObs, & ! IN
894 latLev_T, lonLev_T, & ! IN/OUT
895 latLev_M, lonLev_M ) ! IN/OUT
896
897 latColumn(headerUsedIndex,allkBeg(1)) = latLev_T(1)
898 lonColumn(headerUsedIndex,allkBeg(1)) = lonLev_T(1)
899 end do
900
901 ! gather geographical lat, lon positions of observations from all processors
902 call rpn_comm_allgather(latColumn(:,allkBeg(1)), numHeaderUsedMax, 'MPI_REAL8', &
903 allLatOneLev(:,:), numHeaderUsedMax, 'MPI_REAL8', 'GRID', ierr)
904 call rpn_comm_allgather(lonColumn(:,allkBeg(1)), numHeaderUsedMax, 'MPI_REAL8', &
905 allLonOneLev(:,:), numHeaderUsedMax, 'MPI_REAL8', 'GRID', ierr)
906
907 k_loop: do kIndex = mykBeg, statevector%mykEnd
908 do procIndex = 1, mmpi_nprocs
909 interpInfo%stepProcData(procIndex,stepIndex)%allLat(:,kIndex) = allLatOneLev(1:allNumHeaderUsed(stepIndex,procIndex),procIndex)
910 interpInfo%stepProcData(procIndex,stepIndex)%allLon(:,kIndex) = allLonOneLev(1:allNumHeaderUsed(stepIndex,procIndex),procIndex)
911 end do
912 end do k_loop
913
914 deallocate(latLev_T)
915 deallocate(lonLev_T)
916 deallocate(latLev_M)
917 deallocate(lonLev_M)
918
919 end if ! doSlantPath
920
921 deallocate(lonColumn)
922 deallocate(latColumn)
923
924 end do step_loop2
925
926 if ( gsv_isAllocated(stateVector_1Step) .and. lastCall ) then
927 write(*,*) 's2c_setupInterpInfo: deallocate height3D fields'
928 call gsv_deallocate(stateVector_1Step)
929 end if
930 deallocate(footprintRadiusVec_r4)
931
932 write(*,*) 's2c_setupInterpInfo: latlonChecks and lat/lon MPI comm finished.'
933
934 allocate(allHeaderIndex(numHeaderUsedMax,numStep,mmpi_nprocs))
935 ! gather the headerIndexVec arrays onto all processors
936 call rpn_comm_allgather(headerIndexVec, numHeaderUsedMax*numStep, 'MPI_INTEGER', &
937 allHeaderIndex, numHeaderUsedMax*numStep, 'MPI_INTEGER', &
938 'GRID',ierr)
939
940 do procIndex = 1, mmpi_nprocs
941 do stepIndex = 1, numStep
942 do headerIndex = 1, allNumHeaderUsed(stepIndex,procIndex)
943 interpInfo%stepProcData(procIndex, stepIndex)%allHeaderIndex(headerIndex) = allHeaderIndex(headerIndex,stepIndex,procIndex)
944 end do
945 end do
946 end do
947
948 ! create kdtree to use in footprint operator, if any footprint radius > 0.
949 interpInfo%inputStateVectorType = inputStateVectorType
950 if ( any(allFootprintRadius_r4(:,:,:) > 0.0) ) then
951 if ( (inputStateVectorType == 'nl' .and. .not. associated(tree_nl)) .or. &
952 (inputStateVectorType == 'tl' .and. .not. associated(tree_tlad)) .or. &
953 (inputStateVectorType == 'ad' .and. .not. associated(tree_tlad)) ) then
954
955 write(*,*) 's2c_setupInterpInfo: start creating kdtree for inputStateVectorType=', &
956 inputStateVectorType
957 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
958
959 allocate(positionArray(3,statevector%hco%ni*statevector%hco%nj))
960
961 gridIndex = 0
962 do latIndex = 1, statevector%hco%nj
963 do lonIndex = 1, statevector%hco%ni
964 gridIndex = gridIndex + 1
965 lat = real(stateVector % hco % lat2d_4(lonIndex,latIndex), 8)
966 lon = real(stateVector % hco % lon2d_4(lonIndex,latIndex), 8)
967
968 positionArray(:,gridIndex) = kdtree2_3dPosition(lon, lat)
969
970 end do
971 end do
972
973 nullify(tree)
974 tree => kdtree2_create(positionArray, sort=.false., rearrange=.true.)
975
976 if ( inputStateVectorType == 'nl' ) then
977 tree_nl => tree
978 else
979 tree_tlad => tree
980 end if
981
982 deallocate(positionArray)
983
984 write(*,*) 's2c_setupInterpInfo: done creating kdtree for inputStateVectorType=', &
985 inputStateVectorType
986 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
987
988 end if
989 end if
990
991 do stepIndex = 1, numStep
992 !$OMP PARALLEL DO PRIVATE (procIndex, kIndex, headerIndex, lat_deg_r4, lon_deg_r4, ierr, xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
993 !$OMP subGridIndex, numSubGridsForInterp, subGridForInterp, lat, lon, latRot, lonRot, footprintRadius_r4, numGridpt)
994 do procIndex = 1, mmpi_nprocs
995 do kIndex = mykBeg, statevector%mykEnd
996 do headerIndex = 1, allNumHeaderUsed(stepIndex,procIndex)
997
998 ! Compute the rotated lat/lon, needed for the winds
999
1000 lat_deg_r4 = real(interpInfo%stepProcData(procIndex,stepIndex)%allLat(headerIndex,kIndex) * &
1001 MPC_DEGREES_PER_RADIAN_R8)
1002 lon_deg_r4 = real(interpInfo%stepProcData(procIndex,stepIndex)%allLon(headerIndex,kIndex) * &
1003 MPC_DEGREES_PER_RADIAN_R8)
1004 ierr = gpos_getPositionXY( stateVector%hco%EZscintID, &
1005 xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
1006 lat_deg_r4, lon_deg_r4, subGridIndex )
1007
1008 if ( subGridIndex == 3 ) then
1009 ! both subGrids involved in interpolation, so first treat subGrid 1
1010 numSubGridsForInterp = 2
1011 subGridIndex = 1
1012 else
1013 ! only 1 subGrid involved in interpolation
1014 numSubGridsForInterp = 1
1015 end if
1016
1017 do subGridForInterp = 1, numSubGridsForInterp
1018
1019 if ( subGridForInterp == 1 ) then
1020 ! when only 1 subGrid involved, subGridIndex can be 1 or 2
1021 else
1022 ! when 2 subGrids, subGridIndex is set to 1 for 1st iteration, 2 for second
1023 subGridIndex = 2
1024 end if
1025
1026 if ( interpInfo%hco%rotated .and. &
1027 (gsv_varExist(varName='UU') .or. &
1028 gsv_varExist(varName='VV')) ) then
1029 lat = interpInfo%stepProcData(procIndex,stepIndex)%allLat(headerIndex,kIndex)
1030 lon = interpInfo%stepProcData(procIndex,stepIndex)%allLon(headerIndex,kIndex)
1031 call uvr_RotateLatLon( interpInfo%uvr, & ! INOUT
1032 subGridIndex, & ! IN
1033 latRot, lonRot, & ! OUT (radians)
1034 lat, lon, & ! IN (radians)
1035 'ToLatLonRot') ! IN
1036 interpInfo%stepProcData(procIndex,stepIndex)%allLatRot(subGridIndex,headerIndex,kIndex) = latRot
1037 interpInfo%stepProcData(procIndex,stepIndex)%allLonRot(subGridIndex,headerIndex,kIndex) = lonRot
1038 end if
1039
1040 end do ! subGridForInterp
1041
1042 ! Count total number of grid points for allocating interp depot
1043
1044 footprintRadius_r4 = allFootprintRadius_r4(headerIndex,stepIndex,procIndex)
1045
1046 call s2c_setupHorizInterp(footprintRadius_r4, interpInfo, &
1047 stateVector, headerIndex, kIndex, stepIndex, &
1048 procIndex, numGridpt)
1049
1050 ! for now, just store the number of gridpts for each obs in depotIndexEnd
1051 if ( (subGridIndex == 1) .or. (subGridIndex == 2) ) then
1052 ! indices for only 1 subgrid, other will have zeros
1053 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex,headerIndex,kIndex) = numGridpt(subGridIndex)
1054 else
1055 ! locations on both subGrids will be averaged
1056 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(1,headerIndex,kIndex) = numGridpt(1)
1057 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(2,headerIndex,kIndex) = numGridpt(2)
1058 end if
1059
1060 end do ! headerIndex
1061 end do ! kIndex
1062 end do ! procIndex
1063 !$OMP END PARALLEL DO
1064 end do ! stepIndex
1065
1066 numGridptTotal = 0
1067 do stepIndex = 1, numStep
1068 do procIndex = 1, mmpi_nprocs
1069 do kIndex = mykBeg, statevector%mykEnd
1070 do headerIndex = 1, allNumHeaderUsed(stepIndex,procIndex)
1071 do subGridIndex = 1, interpInfo%hco%numSubGrid
1072 if ( interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex,headerIndex,kIndex) /= -1 ) then
1073 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex,headerIndex,kIndex) = numGridptTotal + 1
1074 numGridptTotal = numGridptTotal + interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex,headerIndex,kIndex)
1075 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex,headerIndex,kIndex) = numGridptTotal
1076 end if
1077 end do ! subGridIndex
1078 end do ! headerIndex
1079 end do ! kIndex
1080 end do ! procIndex
1081 end do ! stepIndex
1082
1083 deallocate(allHeaderIndex)
1084
1085 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1086
1087 ! now that we know the size, allocate main arrays for storing interpolation information
1088 write(*,*) 's2c_setupInterpInfo: numGridptTotal = ', numGridptTotal
1089 allocate( interpInfo%latIndexDepot(numGridptTotal) )
1090 allocate( interpInfo%lonIndexDepot(numGridptTotal) )
1091 allocate( interpInfo%interpWeightDepot(numGridptTotal) )
1092
1093 call utl_tmg_start(33,'------s2c_SetupWeights')
1094 !$OMP PARALLEL DO PRIVATE (procIndex, stepIndex, kIndex, headerIndex, footprintRadius_r4, numGridpt)
1095 do procIndex = 1, mmpi_nprocs
1096 do stepIndex = 1, numStep
1097 do kIndex = mykBeg, statevector%mykEnd
1098
1099 do headerIndex = 1, allNumHeaderUsed(stepIndex,procIndex)
1100
1101 footprintRadius_r4 = allFootprintRadius_r4(headerIndex, stepIndex, procIndex)
1102
1103 call s2c_setupHorizInterp(footprintRadius_r4, interpInfo, stateVector, &
1104 headerIndex, kIndex, stepIndex, procIndex, numGridpt)
1105
1106 end do ! headerIndex
1107
1108 end do ! kIndex
1109 end do ! stepIndex
1110 end do ! procIndex
1111 !$OMP END PARALLEL DO
1112 call utl_tmg_stop(33)
1113
1114 ! reject obs in obsSpaceData if any processor has zero weight
1115 ! called when a mask exists to catch land contaminated ocean obs
1116 if ( stateVector%oceanMask%maskPresent ) then
1117 call s2c_rejectZeroWeightObs(interpInfo,obsSpaceData,mykBeg,stateVector%mykEnd)
1118 end if
1119
1120 ! on the last call, deallocate the tree_nl/tree_tlad
1121 if ( lastCall ) then
1122 if ( inputStateVectorType == 'nl' .and. associated(tree_nl) ) then
1123 call kdtree2_destroy(tree_nl)
1124 else if ( (inputStateVectorType == 'tl' .or. inputStateVectorType == 'ad') .and. &
1125 associated(tree_tlad) ) then
1126 call kdtree2_destroy(tree_tlad)
1127 end if
1128 end if
1129
1130 ! Count the number of TOVS using footprint operator on one level
1131 if ( useFootprintForTovs ) then
1132 numTovsUsingFootprint = 0
1133 numAllTovs = 0
1134 procIndex = mmpi_myid + 1
1135 do stepIndex = 1, numStep
1136 do headerUsedIndex = 1, allNumHeaderUsed(stepIndex,procIndex)
1137 footprintRadius_r4 = allFootprintRadius_r4(headerUsedIndex, stepIndex, procIndex)
1138 headerIndex = headerIndexVec(headerUsedIndex,stepIndex)
1139 codeType = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
1140
1141 if ( tvs_isIdBurpTovs(codeType) ) then
1142 if ( footprintRadius_r4 > 0.0 ) numTovsUsingFootprint = numTovsUsingFootprint + 1
1143 numAllTovs = numAllTovs + 1
1144 end if
1145 end do
1146 end do
1147
1148 if ( numAllTovs > 0 ) then
1149 write(*,'(A,2(I5,A2),F5.1,A)') 's2c_setupInterpInfo: numTovsUsingFootprint/numAllTovs=', &
1150 numTovsUsingFootprint, ' /', numAllTovs, ' (', &
1151 real(numTovsUsingFootprint) / real(numAllTovs) * 100.0, '%)'
1152 end if
1153 end if
1154
1155 deallocate(allFootprintRadius_r4)
1156 deallocate(allLonOneLev)
1157 deallocate(allLatOneLev)
1158
1159 deallocate(headerIndexVec)
1160 deallocate(allNumHeaderUsed)
1161
1162 interpInfo%initialized = .true.
1163
1164 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1165 write(*,*) 's2c_setupInterpInfo: FINISHED'
1166
1167 end subroutine s2c_setupInterpInfo
1168
1169 !---------------------------------------------------------
1170 ! s2c_tl
1171 !---------------------------------------------------------
1172 subroutine s2c_tl(statevector_in, columnAnlInc, columnTrlOnAnlIncLev, obsSpaceData)
1173 !
1174 !:Purpose: Tangent linear version of the horizontal
1175 ! interpolation, used for the increment (or perturbations).
1176 !
1177 implicit none
1178
1179 ! Arguments:
1180 type(struct_gsv), target, intent(in) :: stateVector_in
1181 type(struct_obs) , intent(inout) :: obsSpaceData
1182 type(struct_columnData) , intent(inout) :: columnAnlInc
1183 type(struct_columnData) , intent(inout) :: columnTrlOnAnlIncLev
1184
1185 ! Locals:
1186 type(struct_gsv) :: stateVector_VarsLevs
1187 type(struct_gsv), pointer :: stateVector
1188 integer :: kIndex, kIndex2, levIndex, kCount, stepIndex, numStep, mykEndExtended
1189 integer :: headerIndex, numHeader, numHeaderMax, yourNumHeader
1190 integer :: procIndex, nsize, ierr, headerUsedIndex
1191 real(8) :: weight
1192 real(8), pointer :: allCols_ptr(:,:)
1193 real(pre_incrReal), pointer :: ptr4d(:,:,:,:)
1194 real(pre_incrReal), pointer :: ptr3d_UV(:,:,:)
1195 real(8), allocatable :: cols_hint(:,:,:)
1196 real(8), allocatable :: cols_send(:,:)
1197 real(8), allocatable :: cols_recv(:,:)
1198 real(8), allocatable :: cols_send_1proc(:)
1199 logical :: rejectOutsideObs
1200 character(len=4) :: varName
1201 character(len=4), pointer :: varNames(:)
1202
1203 call utl_tmg_start(30,'--StateToColumn')
1204
1205 if ( mmpi_myid == 0 ) write(*,*) 's2c_tl: Horizontal interpolation StateVector --> ColumnData'
1206 call utl_tmg_start(38,'----s2c_TL')
1207
1208 call rpn_comm_barrier('GRID',ierr)
1209
1210 if ( .not. gsv_isAllocated(stateVector_in) ) then
1211 call utl_abort('s2c_tl: stateVector must be allocated')
1212 end if
1213
1214 if (interpInfo_tlad%initialized) then
1215 if (.not. hco_equal(interpInfo_tlad%hco,stateVector_in%hco)) then
1216 write(*,*) 's2c_tl: WARNING! Current hco grid parameters differ from allocated interpInfo_tlad!'
1217 write(*,*) 's2c_tl: InterpInfo_tlad will be deallocated.'
1218 call s2c_deallocInterpInfo(inputStateVectorType='tlad')
1219 end if
1220 end if
1221
1222 ! check the column and statevector have same nk/varNameList
1223 call checkColumnStatevectorMatch(columnAnlInc,statevector_in)
1224
1225 ! if we only compute Height and Pressure on column, make copy without them
1226 if ( calcHeightPressIncrOnColumn ) then
1227 allocate(stateVector)
1228 call gsv_allocate( stateVector, statevector_in%numstep, &
1229 statevector_in%hco, statevector_in%vco, &
1230 mpi_local_opt=.true., &
1231 dataKind_opt=gsv_getDataKind(statevector_in), &
1232 allocHeight_opt=.false., allocPressure_opt=.false. )
1233 call gsv_copy(stateVector_in, stateVector, allowVarMismatch_opt=.true.)
1234 else
1235 stateVector => stateVector_in
1236
1237 ! calculate delP_T/delP_M and del Z_T/Z_M on the grid
1238 call gvt_transform( statevector, 'ZandP_tl' )
1239 end if
1240
1241 nullify(varNames)
1242 call gsv_varNamesList(varNames, statevector)
1243 if (statevector%mpi_distribution == 'None') then
1244 call gsv_allocate( statevector_VarsLevs, statevector%numstep, &
1245 statevector%hco, statevector%vco, &
1246 mpi_local_opt=.false., mpi_distribution_opt='None', &
1247 dataKind_opt=gsv_getDataKind(statevector), &
1248 varNames_opt=varNames )
1249 else
1250 call gsv_allocate( statevector_VarsLevs, statevector%numstep, &
1251 statevector%hco, statevector%vco, &
1252 mpi_local_opt=.true., mpi_distribution_opt='VarsLevs', &
1253 dataKind_opt=gsv_getDataKind(statevector), &
1254 varNames_opt=varNames )
1255 end if
1256 deallocate(varNames)
1257 call gsv_transposeTilesToVarsLevs( statevector, statevector_VarsLevs )
1258
1259 numStep = stateVector_VarsLevs%numStep
1260 numHeader = obs_numheader(obsSpaceData)
1261 call rpn_comm_allreduce(numHeader, numHeaderMax, 1, &
1262 'MPI_INTEGER', 'MPI_MAX', 'GRID', ierr)
1263
1264 if ( .not. interpInfo_tlad%initialized ) then
1265 rejectOutsideObs = .false.
1266 call utl_tmg_stop(38)
1267 call utl_tmg_start(31,'----s2c_Setups')
1268 call s2c_setupInterpInfo( interpInfo_tlad, obsSpaceData, stateVector_VarsLevs, &
1269 1, numHeader, timeInterpType_tlad, rejectOutsideObs, &
1270 inputStateVectorType='tl' )
1271 call utl_tmg_stop(31)
1272 call utl_tmg_start(38,'----s2c_TL')
1273 end if
1274
1275 ! arrays for interpolated column for 1 level/variable and each time step
1276 allocate(cols_hint(maxval(interpInfo_tlad%allNumHeaderUsed),numStep,mmpi_nprocs))
1277 cols_hint(:,:,:) = 0.0d0
1278
1279 ! arrays for sending/receiving time interpolated column for 1 level/variable
1280 allocate(cols_send(numHeaderMax,mmpi_nprocs))
1281 cols_send(:,:) = 0.0d0
1282
1283 allocate(cols_recv(numHeaderMax,mmpi_nprocs))
1284 cols_recv(:,:) = 0.0d0
1285
1286 allocate(cols_send_1proc(numHeaderMax))
1287
1288 ! set contents of column to zero
1289 allCols_ptr => col_getAllColumns(columnAnlInc)
1290 if ( numHeader > 0 ) allCols_ptr(:,:) = 0.0d0
1291
1292 call gsv_getField(stateVector_VarsLevs, ptr4d)
1293
1294 mykEndExtended = stateVector_VarsLevs%mykBeg + maxval(stateVector_VarsLevs%allkCount(:)) - 1
1295
1296 kCount = 0
1297 k_loop: do kIndex = stateVector_VarsLevs%mykBeg, mykEndExtended
1298
1299 kCount = kCount + 1
1300
1301 if ( kIndex <= stateVector_VarsLevs%mykEnd ) then
1302 varName = gsv_getVarNameFromK(statevector,kIndex)
1303
1304 if ( varName == 'UU' .or. varName == 'VV' ) then
1305 call gsv_getFieldUV(stateVector_VarsLevs,ptr3d_UV,kIndex)
1306 end if
1307
1308 call utl_tmg_start(39,'------s2c_TL_Hinterp')
1309 !$OMP PARALLEL DO PRIVATE (stepIndex, procIndex, yourNumHeader, headerIndex)
1310 step_loop: do stepIndex = 1, numStep
1311 if ( maxval(interpInfo_tlad%allNumHeaderUsed(stepIndex,:)) == 0 ) cycle step_loop
1312
1313 ! interpolate to the columns destined for all procs for all steps and one lev/var
1314 do procIndex = 1, mmpi_nprocs
1315 yourNumHeader = interpInfo_tlad%allNumHeaderUsed(stepIndex,procIndex)
1316 if ( yourNumHeader > 0 ) then
1317 if ( varName == 'UU' ) then
1318 call myezuvint_tl( cols_hint(1:yourNumHeader,stepIndex,procIndex), 'UU', &
1319 ptr4d(:,:,kIndex,stepIndex), ptr3d_UV(:,:,stepIndex), &
1320 interpInfo_tlad, kIndex, stepIndex, procIndex )
1321 else if ( varName == 'VV' ) then
1322 call myezuvint_tl( cols_hint(1:yourNumHeader,stepIndex,procIndex), 'VV', &
1323 ptr3d_UV(:,:,stepIndex), ptr4d(:,:,kIndex,stepIndex), &
1324 interpInfo_tlad, kIndex, stepIndex, procIndex )
1325 else
1326 call myezsint_tl( cols_hint(1:yourNumHeader,stepIndex,procIndex), &
1327 ptr4d(:,:,kIndex,stepIndex), interpInfo_tlad, kIndex, &
1328 stepIndex, procIndex )
1329 end if
1330 end if
1331 end do
1332
1333 end do step_loop
1334 !$OMP END PARALLEL DO
1335 call utl_tmg_stop(39)
1336
1337 ! interpolate in time to the columns destined for all procs and one level/variable
1338 do procIndex = 1, mmpi_nprocs
1339 cols_send_1proc(:) = 0.0d0
1340 do stepIndex = 1, numStep
1341 !$OMP PARALLEL DO PRIVATE (headerUsedIndex, headerIndex, weight)
1342 do headerUsedIndex = 1, interpInfo_tlad%allNumHeaderUsed(stepIndex, procIndex)
1343 headerIndex = interpInfo_tlad%stepProcData(procIndex, stepIndex)%allHeaderIndex(headerUsedIndex)
1344 weight = oti_getTimeInterpWeightMpiGlobal(interpInfo_tlad%oti, &
1345 headerIndex,stepIndex,procIndex)
1346 cols_send_1proc(headerIndex) = cols_send_1proc(headerIndex) &
1347 + weight * cols_hint(headerUsedIndex,stepIndex,procIndex)
1348
1349 end do
1350 !$OMP END PARALLEL DO
1351 end do
1352 cols_send(:,procIndex) = cols_send_1proc(:)
1353 end do
1354
1355 else
1356
1357 ! this value of k does not exist on this mpi task
1358 cols_send(:,:) = 0.0
1359
1360 end if ! if kIndex <= mykEnd
1361
1362 call rpn_comm_barrier('GRID',ierr)
1363
1364 ! mpi communication: alltoall for one level/variable
1365 nsize = numHeaderMax
1366 if(mmpi_nprocs > 1) then
1367 call rpn_comm_alltoall(cols_send, nsize, 'MPI_REAL8', &
1368 cols_recv, nsize, 'MPI_REAL8', 'GRID', ierr)
1369 else
1370 cols_recv(:,1) = cols_send(:,1)
1371 end if
1372
1373 ! reorganize ensemble of distributed columns
1374 !$OMP PARALLEL DO PRIVATE (procIndex, kIndex2, varName, levIndex, allCols_ptr, headerIndex)
1375 proc_loop: do procIndex = 1, mmpi_nprocs
1376 ! This is kIndex value of source (can be different for destination)
1377 kIndex2 = statevector_VarsLevs%allkBeg(procIndex) + kCount - 1
1378 if ( kIndex2 > stateVector_VarsLevs%allkEnd(procIndex) ) cycle proc_loop
1379
1380 ! Figure out which variable/level of destination
1381 varName = gsv_getVarNameFromK(statevector,kIndex2)
1382 levIndex = gsv_getLevFromK(statevector,kIndex2)
1383 allCols_ptr => col_getAllColumns(columnAnlInc,varName)
1384
1385 do headerIndex = 1, numHeader
1386 allCols_ptr(levIndex,headerIndex) = cols_recv(headerIndex,procIndex)
1387 end do
1388 end do proc_loop
1389 !$OMP END PARALLEL DO
1390
1391 end do k_loop
1392
1393 if (calcHeightPressIncrOnColumn) then
1394 ! calculate delP_T/delP_M and del Z_T/Z_M on the columns
1395 call czp_calcZandP_tl(columnAnlInc, columnTrlOnAnlIncLev)
1396 end if
1397
1398 deallocate(cols_hint)
1399 deallocate(cols_send)
1400 deallocate(cols_recv)
1401 deallocate(cols_send_1proc)
1402
1403 call gsv_deallocate( statevector_VarsLevs )
1404 if (calcHeightPressIncrOnColumn) then
1405 call gsv_deallocate( stateVector )
1406 deallocate(stateVector)
1407 end if
1408
1409 if (slantPath_TO_tlad) call pressureProfileMonotonicityCheck(obsSpaceData, columnTrlOnAnlIncLev)
1410
1411 call utl_tmg_stop(38)
1412
1413 call utl_tmg_stop(30)
1414
1415 end subroutine s2c_tl
1416
1417 !---------------------------------------------------------
1418 ! s2c_ad
1419 !---------------------------------------------------------
1420 subroutine s2c_ad(statevector_out, columnAnlInc, columnTrlOnAnlIncLev, obsSpaceData)
1421 !
1422 !:Purpose: Adjoint version of the horizontal interpolation,
1423 ! used for the cost function gradient with respect to the increment.
1424 !
1425 implicit none
1426
1427 ! Arguments:
1428 type(struct_gsv), target, intent(inout) :: stateVector_out
1429 type(struct_obs) , intent(inout) :: obsSpaceData
1430 type(struct_columnData) , intent(inout) :: columnAnlInc
1431 type(struct_columnData) , intent(inout) :: columnTrlOnAnlIncLev
1432
1433 ! Locals:
1434 type(struct_gsv) :: stateVector_VarsLevs
1435 type(struct_gsv), pointer :: stateVector
1436 integer :: kIndex, kIndex2, kCount, levIndex, stepIndex, numStep, mykEndExtended
1437 integer :: headerIndex, numHeader, numHeaderMax, yourNumHeader
1438 integer :: procIndex, nsize, ierr, headerUsedIndex
1439 character(len=4) :: varName
1440 real(8) :: weight
1441 real(8), pointer :: allCols_ptr(:,:)
1442 real(pre_incrReal), pointer :: ptr4d(:,:,:,:), ptr3d_UV(:,:,:)
1443 real(8), allocatable :: cols_hint(:,:,:)
1444 real(8), allocatable :: cols_send(:,:)
1445 real(8), allocatable :: cols_recv(:,:)
1446 logical :: rejectOutsideObs
1447 character(len=4), pointer :: varNames(:)
1448
1449 call utl_tmg_start(30,'--StateToColumn')
1450
1451 if(mmpi_myid == 0) write(*,*) 's2c_ad: Adjoint of horizontal interpolation StateVector --> ColumnData'
1452 call utl_tmg_start(40,'----s2c_AD')
1453
1454 call rpn_comm_barrier('GRID',ierr)
1455
1456 if ( .not. gsv_isAllocated(stateVector_out) ) then
1457 call utl_abort('s2c_ad: stateVector must be allocated')
1458 end if
1459
1460 if (interpInfo_tlad%initialized) then
1461 if (.not. hco_equal(interpInfo_tlad%hco,stateVector_out%hco)) then
1462 write(*,*) 's2c_ad: WARNING! Current hco grid parameters differ from allocated interpInfo_tlad!'
1463 write(*,*) 's2c_ad: InterpInfo_tlad will be deallocated.'
1464 call s2c_deallocInterpInfo(inputStateVectorType='tlad')
1465 end if
1466 end if
1467
1468
1469 ! if we only compute Height and Pressure on column, make copy without them
1470 if (calcHeightPressIncrOnColumn) then
1471 allocate(stateVector)
1472 call gsv_allocate( stateVector, statevector_out%numstep, &
1473 statevector_out%hco, statevector_out%vco, &
1474 mpi_local_opt=.true., &
1475 dataKind_opt=gsv_getDataKind(statevector_out), &
1476 allocHeight_opt=.false., allocPressure_opt=.false. )
1477 ! Adjoint of calculate del Z_T/Z_M and delP_T/delP_M on the columns
1478 call czp_calcZandP_ad(columnAnlInc, columnTrlOnAnlIncLev)
1479 else
1480 stateVector => stateVector_out
1481 end if
1482
1483 nullify(varNames)
1484 call gsv_varNamesList(varNames, statevector)
1485 call gsv_allocate( statevector_VarsLevs, statevector%numstep, &
1486 statevector%hco, statevector%vco, &
1487 mpi_local_opt=.true., mpi_distribution_opt='VarsLevs', &
1488 dataKind_opt=gsv_getDataKind(statevector), &
1489 varNames_opt=varNames )
1490 deallocate(varNames)
1491 call gsv_zero( statevector_VarsLevs )
1492
1493 numStep = stateVector_VarsLevs%numStep
1494 numHeader = obs_numheader(obsSpaceData)
1495 call rpn_comm_allreduce(numHeader, numHeaderMax, 1, &
1496 'MPI_INTEGER', 'MPI_MAX', 'GRID', ierr)
1497
1498 if ( .not. interpInfo_tlad%initialized ) then
1499 rejectOutsideObs = .false.
1500 call utl_tmg_stop(40)
1501 call utl_tmg_start(31,'----s2c_Setups')
1502 call s2c_setupInterpInfo( interpInfo_tlad, obsSpaceData, stateVector_VarsLevs, &
1503 1, numHeader, timeInterpType_tlad, rejectOutsideObs, &
1504 inputStateVectorType='ad' )
1505 call utl_tmg_stop(31)
1506 call utl_tmg_start(40,'----s2c_AD')
1507 end if
1508
1509 ! arrays for interpolated column for 1 level/variable and each time step
1510 allocate(cols_hint(maxval(interpInfo_tlad%allNumHeaderUsed),numStep,mmpi_nprocs))
1511 cols_hint(:,:,:) = 0.0d0
1512
1513 ! arrays for sending/receiving time interpolated column for 1 level/variable
1514 allocate(cols_send(numHeaderMax,mmpi_nprocs))
1515 cols_send(:,:) = 0.0d0
1516
1517 allocate(cols_recv(numHeaderMax,mmpi_nprocs))
1518 cols_recv(:,:) = 0.0d0
1519
1520 ! set contents of column to zero
1521 allCols_ptr => col_getAllColumns(columnAnlInc)
1522
1523 call gsv_getField(stateVector_VarsLevs,ptr4d)
1524 mykEndExtended = stateVector_VarsLevs%mykBeg + maxval(stateVector_VarsLevs%allkCount(:)) - 1
1525
1526 kCount = 0
1527 k_loop: do kIndex = stateVector_VarsLevs%mykBeg, mykEndExtended
1528
1529 kCount = kCount + 1
1530
1531 ! reorganize ensemble of distributed columns
1532 !$OMP PARALLEL DO PRIVATE (procIndex, kIndex2, varName, levIndex, allCols_ptr, headerIndex)
1533 proc_loop: do procIndex = 1, mmpi_nprocs
1534 ! This is kIndex value of destination (can be different for source)
1535 kIndex2 = statevector_VarsLevs%allkBeg(procIndex) + kCount - 1
1536 if ( kIndex2 > stateVector_VarsLevs%allkEnd(procIndex) ) cycle proc_loop
1537
1538 ! Figure out which variable/level of source
1539 varName = gsv_getVarNameFromK(statevector,kIndex2)
1540 levIndex = gsv_getLevFromK(statevector,kIndex2)
1541 allCols_ptr => col_getAllColumns(columnAnlInc,varName)
1542
1543 do headerIndex = 1, numHeader
1544 cols_send(headerIndex,procIndex) = allCols_ptr(levIndex,headerIndex)
1545 end do
1546 end do proc_loop
1547 !$OMP END PARALLEL DO
1548
1549 call rpn_comm_barrier('GRID',ierr)
1550
1551 ! mpi communication: alltoall for one level/variable
1552 nsize = numHeaderMax
1553 if(mmpi_nprocs > 1) then
1554 call rpn_comm_alltoall(cols_send, nsize, 'MPI_REAL8', &
1555 cols_recv, nsize, 'MPI_REAL8', 'GRID', ierr)
1556 else
1557 cols_recv(:,1) = cols_send(:,1)
1558 end if
1559
1560 if ( kIndex <= stateVector_VarsLevs%mykEnd ) then
1561 varName = gsv_getVarNameFromK(statevector,kIndex)
1562
1563 if ( varName == 'UU' .or. varName == 'VV' ) then
1564 call gsv_getFieldUV(stateVector_VarsLevs, ptr3d_UV, kIndex)
1565 end if
1566
1567 ! interpolate in time to the columns destined for all procs and one level/variable
1568 do procIndex = 1, mmpi_nprocs
1569 do stepIndex = 1, numStep
1570 !$OMP PARALLEL DO PRIVATE (headerIndex, headerUsedIndex, weight)
1571 do headerUsedIndex = 1, interpInfo_tlad%allNumHeaderUsed(stepIndex, procIndex)
1572
1573 headerIndex = interpInfo_tlad%stepProcData(procIndex, stepIndex)%allHeaderIndex(headerUsedIndex)
1574 weight = oti_getTimeInterpWeightMpiGlobal(interpInfo_tlad%oti, &
1575 headerIndex,stepIndex,procIndex)
1576
1577 cols_hint(headerUsedIndex,stepIndex,procIndex) = &
1578 weight * cols_recv(headerIndex,procIndex)
1579
1580 end do
1581 !$OMP END PARALLEL DO
1582 end do
1583 end do
1584
1585 call utl_tmg_start(41,'------s2c_AD_Hinterp')
1586 !$OMP PARALLEL DO PRIVATE (stepIndex, procIndex, yourNumHeader)
1587 step_loop: do stepIndex = 1, numStep
1588 if ( maxval(interpInfo_tlad%allNumHeaderUsed(stepIndex,:)) == 0 ) cycle step_loop
1589
1590 ! interpolate to the columns destined for all procs for all steps and one lev/var
1591 do procIndex = 1, mmpi_nprocs
1592 yourNumHeader = interpInfo_tlad%allNumHeaderUsed(stepIndex,procIndex)
1593 if ( yourNumHeader > 0 ) then
1594 if ( varName == 'UU' ) then
1595 call myezuvint_ad( cols_hint(1:yourNumHeader,stepIndex,procIndex), 'UU', &
1596 ptr4d(:,:,kIndex,stepIndex), ptr3d_UV(:,:,stepIndex), &
1597 interpInfo_tlad, kIndex, stepIndex, procIndex )
1598 else if ( varName == 'VV' ) then
1599 call myezuvint_ad( cols_hint(1:yourNumHeader,stepIndex,procIndex), 'VV', &
1600 ptr3d_UV(:,:,stepIndex), ptr4d(:,:,kIndex,stepIndex), &
1601 interpInfo_tlad, kIndex, stepIndex, procIndex )
1602 else
1603 call myezsint_ad( cols_hint(1:yourNumHeader,stepIndex,procIndex), &
1604 ptr4d(:,:,kIndex,stepIndex), interpInfo_tlad, kIndex, &
1605 stepIndex, procIndex )
1606 end if
1607 end if
1608 end do
1609
1610 end do step_loop
1611 !$OMP END PARALLEL DO
1612 call utl_tmg_stop(41)
1613
1614 end if ! if kIndex <= mykEnd
1615
1616 end do k_loop
1617
1618 deallocate(cols_hint)
1619 deallocate(cols_send)
1620 deallocate(cols_recv)
1621
1622 call rpn_comm_barrier('GRID',ierr)
1623
1624 call gsv_transposeTilesToVarsLevsAd( statevector_VarsLevs, statevector )
1625
1626 if (calcHeightPressIncrOnColumn) then
1627 call gsv_zero(statevector_out)
1628 call gsv_copy(stateVector, stateVector_out, allowVarMismatch_opt=.true.)
1629 else
1630 ! Adjoint of calculate del Z_T/Z_M and delP_T/delP_M on the grid
1631 call gvt_transform( statevector, 'ZandP_ad' )
1632 end if
1633
1634 call gsv_deallocate( statevector_VarsLevs )
1635 if (calcHeightPressIncrOnColumn) then
1636 call gsv_deallocate( statevector )
1637 deallocate(stateVector)
1638 end if
1639
1640 if (slantPath_TO_tlad) call pressureProfileMonotonicityCheck(obsSpaceData, columnTrlOnAnlIncLev)
1641
1642 call utl_tmg_stop(40)
1643
1644 call utl_tmg_stop(30)
1645
1646 end subroutine s2c_ad
1647
1648 !---------------------------------------------------------
1649 ! s2c_nl
1650 !---------------------------------------------------------
1651 subroutine s2c_nl(stateVector, obsSpaceData, column, hco_core, timeInterpType, &
1652 varName_opt, numObsBatches_opt, dealloc_opt, moveObsAtPole_opt, &
1653 beSilent_opt)
1654 !
1655 !:Purpose: Non-linear version of the horizontal interpolation,
1656 ! used for a full field (usually the background state when computing
1657 ! the innovation vector).
1658 !
1659 implicit none
1660
1661 ! Arguments:
1662 type(struct_gsv) , intent(inout) :: stateVector
1663 type(struct_obs) , intent(inout) :: obsSpaceData
1664 type(struct_columnData), intent(inout) :: column
1665 type(struct_hco), pointer, intent(in) :: hco_core
1666 character(len=*) , intent(in) :: timeInterpType
1667 character(len=*), optional, intent(in) :: varName_opt
1668 integer, optional, intent(in) :: numObsBatches_opt
1669 logical, optional, intent(in) :: dealloc_opt
1670 logical, optional, intent(in) :: moveObsAtPole_opt
1671 logical, optional, intent(in) :: beSilent_opt
1672
1673 ! Locals:
1674 type(struct_gsv), save :: stateVector_VarsLevs
1675 integer :: kIndex, kIndex2, kCount, stepIndex, numStep, mykEndExtended
1676 integer :: headerIndex, headerIndex2, numHeader, numHeaderMax, yourNumHeader
1677 integer :: headerIndexBeg, headerIndexEnd, obsBatchIndex, numObsBatches
1678 integer :: procIndex, nsize, ierr, headerUsedIndex, allHeaderIndexBeg(mmpi_nprocs)
1679 integer :: kIndexHeightSfc, varNameIndex, allNumHeader(mmpi_nprocs)
1680 real(8) :: weight
1681 character(len=4) :: varName
1682 real(8), pointer :: column_ptr(:), ptr2d_r8(:,:), allCols_ptr(:,:)
1683 real(4), pointer :: ptr4d_r4(:,:,:,:), ptr3d_UV_r4(:,:,:)
1684 real(8), allocatable :: cols_hint(:,:,:)
1685 real(8), allocatable :: cols_send(:,:)
1686 real(8), allocatable :: cols_recv(:,:)
1687 real(8), allocatable :: cols_send_1proc(:)
1688 integer :: displs(mmpi_nprocs), nsizes(mmpi_nprocs)
1689 integer :: senddispls(mmpi_nprocs), sendsizes(mmpi_nprocs)
1690 integer :: recvdispls(mmpi_nprocs), recvsizes(mmpi_nprocs)
1691 logical :: dealloc, moveObsAtPole, rejectOutsideObs, beSilent
1692 logical, save :: firstCall = .true.
1693 character(len=4), pointer :: varNames(:)
1694
1695 call utl_tmg_start(30,'--StateToColumn')
1696 call utl_tmg_start(34,'----s2c_NL')
1697
1698 call utl_tmg_start(37,'------s2c_NL_barrier')
1699 call rpn_comm_barrier('GRID',ierr)
1700 call utl_tmg_stop(37)
1701
1702 if ( present(beSilent_opt) ) then
1703 beSilent = beSilent_opt
1704 else
1705 beSilent = .false.
1706 end if
1707
1708 if ( .not. beSilent ) then
1709 write(*,*) 's2c_nl: STARTING'
1710 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1711 end if
1712
1713 if ( .not. gsv_isAllocated(stateVector) ) then
1714 call utl_abort('s2c_nl: stateVector must be allocated')
1715 end if
1716
1717 if (present(dealloc_opt)) then
1718 dealloc = dealloc_opt
1719 else
1720 dealloc = .true.
1721 end if
1722
1723 ! determine number of obs batches (to reduce memory usage)
1724 if (present(numObsBatches_opt)) then
1725 numObsBatches = numObsBatches_opt
1726 else
1727 numObsBatches = 1
1728 end if
1729 if (.not. dealloc) then
1730 numObsBatches = 1 ! multiple batches only possible if dealloc=.true.
1731 end if
1732
1733 if (interpInfo_nl%initialized) then
1734 if (.not. hco_equal(interpInfo_nl%hco,stateVector%hco) .or. numObsBatches > 1) then
1735 write(*,*) 's2c_nl: WARNING! Current hco grid parameters differ from allocated interpInfo!'
1736 write(*,*) 's2c_nl: InterpInfo will be deallocated.'
1737 call s2c_deallocInterpInfo(inputStateVectorType='nl')
1738 end if
1739 end if
1740
1741 if ( stateVector%mpi_distribution /= 'Tiles' ) then
1742 call utl_abort('s2c_nl: stateVector must by Tiles distributed')
1743 end if
1744
1745 if ( present(moveObsAtPole_opt) ) then
1746 moveObsAtPole = moveObsAtPole_opt
1747 else
1748 moveObsAtPole = .false.
1749 end if
1750
1751 ! check the column and statevector have same nk/varNameList
1752 call checkColumnStatevectorMatch(column,statevector)
1753
1754 ! calculate delP_T/delP_M and del Z_T/Z_M on the grid
1755 call gvt_transform( statevector, 'ZandP_nl' )
1756
1757 if ( dealloc .or. firstCall ) then
1758 nullify(varNames)
1759 call gsv_varNamesList(varNames, statevector)
1760 call gsv_allocate( statevector_VarsLevs, stateVector%numstep, &
1761 stateVector%hco, stateVector%vco, mpi_local_opt=.true., &
1762 mpi_distribution_opt='VarsLevs', dataKind_opt=4, &
1763 allocHeightSfc_opt=.true., varNames_opt=varNames )
1764 deallocate(varNames)
1765 else
1766 if (mmpi_myid == 0 .and. .not. beSilent) write(*,*) 's2c_nl: avoid re-allocating statevector_VarsLevs'
1767 call gsv_zero(statevector_VarsLevs)
1768 end if
1769
1770 call gsv_transposeTilesToVarsLevs( stateVector, stateVector_VarsLevs, &
1771 beSilent_opt=beSilent )
1772
1773 numStep = stateVector_VarsLevs%numStep
1774
1775 if ( .not. interpInfo_nl%initialized ) then
1776 call utl_tmg_stop(34)
1777 call utl_tmg_start(31,'----s2c_Setups')
1778 ! also reject obs outside (LAM) domain and optionally move obs near
1779 ! numerical pole to first/last analysis grid latitude
1780 call latlonChecksAnlGrid( obsSpaceData, hco_core, moveObsAtPole )
1781
1782 ! Do not reject obs for global domain
1783 rejectOutsideObs = .not. stateVector_VarsLevs%hco%global
1784 write(*,*) 's2c_nl: rejectOutsideObs = ', rejectOutsideObs
1785 call utl_tmg_stop(31)
1786 call utl_tmg_start(34,'----s2c_NL')
1787
1788 end if
1789
1790 ! set contents of column to zero (1 variable or all)
1791 allCols_ptr => col_getAllColumns(column,varName_opt)
1792 if ( obs_numHeader(obsSpaceData) > 0 ) allCols_ptr(:,:) = 0.0d0
1793
1794 OBSBATCH: do obsBatchIndex = 1, numObsBatches
1795 headerIndexBeg = 1 + (obsBatchIndex - 1) * (obs_numheader(obsSpaceData) / numObsBatches)
1796 if (obsBatchIndex == numObsBatches) then
1797 headerIndexEnd = obs_numheader(obsSpaceData)
1798 else
1799 headerIndexEnd = headerIndexBeg + (obs_numheader(obsSpaceData) / numObsBatches) - 1
1800 end if
1801 numHeader = headerIndexEnd - headerIndexBeg + 1
1802 call rpn_comm_allreduce(numHeader, numHeaderMax, 1, &
1803 'MPI_INTEGER', 'MPI_MAX', 'GRID', ierr)
1804 if ( .not. beSilent ) then
1805 write(*,*) 's2c_nl: headerIndexBeg/End, numHeader, numHeaderMax = ', &
1806 headerIndexBeg, headerIndexEnd, numHeader, numHeaderMax
1807 if (mmpi_myid == 0) then
1808 write(*,*) 's2c_nl: min/max of allNumHeader = ', &
1809 minval(allNumHeader), maxval(allNumHeader)
1810 end if
1811 end if
1812
1813 call rpn_comm_allgather(numHeader, 1,'mpi_integer', &
1814 allNumHeader,1,'mpi_integer','grid',ierr)
1815 call rpn_comm_allgather(headerIndexBeg, 1,'mpi_integer', &
1816 allHeaderIndexBeg,1,'mpi_integer','grid',ierr)
1817
1818 if ( .not. interpInfo_nl%initialized ) then
1819 call utl_tmg_stop(34)
1820 call utl_tmg_start(31,'----s2c_Setups')
1821
1822 ! compute and collect all obs grids onto all mpi tasks
1823 call s2c_setupInterpInfo( interpInfo_nl, obsSpaceData, stateVector_VarsLevs, &
1824 headerIndexBeg, headerIndexEnd, &
1825 timeInterpType, rejectOutsideObs, &
1826 inputStateVectorType='nl', &
1827 lastCall_opt=(obsBatchIndex==numObsBatches))
1828 if ( mmpi_myid == 0 .and. verbose ) then
1829 do stepIndex = 1, numStep
1830 write(*,*) 's2c_nl: stepIndex, allNumHeaderUsed = ', &
1831 stepIndex, interpInfo_nl%allNumHeaderUsed(stepIndex,:)
1832 end do
1833 end if
1834 call utl_tmg_stop(31)
1835 call utl_tmg_start(34,'----s2c_NL')
1836 end if
1837
1838 ! arrays for interpolated column for 1 level/variable and each time step
1839 allocate(cols_hint(maxval(interpInfo_nl%allNumHeaderUsed),numStep,mmpi_nprocs))
1840 cols_hint(:,:,:) = 0.0d0
1841
1842 ! arrays for sending/receiving time interpolated column for 1 level/variable
1843 allocate(cols_send(numHeaderMax,mmpi_nprocs))
1844 cols_send(:,:) = 0.0d0
1845
1846 allocate(cols_recv(numHeaderMax,mmpi_nprocs))
1847 cols_recv(:,:) = 0.0d0
1848
1849 allocate(cols_send_1proc(numHeaderMax))
1850
1851 call gsv_getField(stateVector_VarsLevs,ptr4d_r4)
1852
1853 mykEndExtended = stateVector_VarsLevs%mykBeg + maxval(stateVector_VarsLevs%allkCount(:)) - 1
1854
1855 kCount = 0
1856 k_loop: do kIndex = stateVector_VarsLevs%mykBeg, mykEndExtended
1857 kCount = kCount + 1
1858
1859 if ( kIndex <= stateVector_VarsLevs%mykEnd ) then
1860 varName = gsv_getVarNameFromK(stateVector_VarsLevs,kIndex)
1861
1862 call utl_tmg_start(35,'------s2c_NL_Hinterp')
1863 if ( varName == 'UU' .or. varName == 'VV' ) then
1864 call gsv_getFieldUV(stateVector_VarsLevs,ptr3d_UV_r4,kIndex)
1865 end if
1866
1867 step_loop: do stepIndex = 1, numStep
1868 if ( maxval(interpInfo_nl%allNumHeaderUsed(stepIndex,:)) == 0 ) cycle step_loop
1869
1870 ! interpolate to the columns destined for all procs for all steps and one lev/var
1871 !$OMP PARALLEL DO PRIVATE (procIndex, yourNumHeader)
1872 do procIndex = 1, mmpi_nprocs
1873 yourNumHeader = interpInfo_nl%allNumHeaderUsed(stepIndex,procIndex)
1874 if ( yourNumHeader > 0 ) then
1875 if ( varName == 'UU' ) then
1876 call myezuvint_nl( cols_hint(1:yourNumHeader,stepIndex,procIndex), 'UU', &
1877 ptr4d_r4(:,:,kIndex,stepIndex), ptr3d_UV_r4(:,:,stepIndex), &
1878 interpInfo_nl, kindex, stepIndex, procIndex )
1879 else if ( varName == 'VV' ) then
1880 call myezuvint_nl( cols_hint(1:yourNumHeader,stepIndex,procIndex), 'VV', &
1881 ptr3d_UV_r4(:,:,stepIndex), ptr4d_r4(:,:,kIndex,stepIndex), &
1882 interpInfo_nl, kindex, stepIndex, procIndex )
1883 else
1884 call myezsint_nl( cols_hint(1:yourNumHeader,stepIndex,procIndex), &
1885 ptr4d_r4(:,:,kIndex,stepIndex), &
1886 interpInfo_nl, kindex, stepIndex, procIndex )
1887 end if
1888 end if
1889 end do
1890 !$OMP END PARALLEL DO
1891
1892 end do step_loop
1893 call utl_tmg_stop(35)
1894
1895 ! interpolate in time to the columns destined for all procs and one level/variable
1896 do procIndex = 1, mmpi_nprocs
1897 cols_send_1proc(:) = 0.0d0
1898 do stepIndex = 1, numStep
1899 !$OMP PARALLEL DO PRIVATE (headerIndex, headerIndex2, headerUsedIndex, weight)
1900 do headerUsedIndex = 1, interpInfo_nl%allNumHeaderUsed(stepIndex, procIndex)
1901 headerIndex = interpInfo_nl%stepProcData(procIndex,stepIndex)%allHeaderIndex(headerUsedIndex)
1902 headerIndex2 = headerIndex - allHeaderIndexBeg(procIndex) + 1
1903 weight = oti_getTimeInterpWeightMpiGlobal(interpInfo_nl%oti, &
1904 headerIndex2,stepIndex,procIndex)
1905 cols_send_1proc(headerIndex2) = cols_send_1proc(headerIndex2) + &
1906 weight * cols_hint(headerUsedIndex,stepIndex,procIndex)
1907 end do
1908 !$OMP END PARALLEL DO
1909 end do
1910 cols_send(:,procIndex) = cols_send_1proc(:)
1911 end do
1912
1913 else
1914
1915 ! this value of k does not exist on this mpi task
1916 cols_send(:,:) = 0.0d0
1917
1918 end if ! if kIndex <= mykEnd
1919
1920 call utl_tmg_start(37,'------s2c_NL_barrier')
1921 call rpn_comm_barrier('GRID',ierr)
1922 call utl_tmg_stop(37)
1923
1924 ! mpi communication: alltoallv for one level/variable
1925 call utl_tmg_start(36,'------s2c_NL_allToAll')
1926
1927 ! only receive the data from tasks with data, same amount from all of those
1928 recvsizes(:) = 0
1929 do procIndex = 1, mmpi_nprocs
1930 kIndex2 = stateVector_VarsLevs%allkBeg(procIndex) + kCount - 1
1931 if (kIndex2 > stateVector_VarsLevs%allkEnd(procIndex)) cycle
1932 recvsizes(procIndex) = numHeader
1933 end do
1934 recvdispls(1) = 0
1935 do procIndex = 2, mmpi_nprocs
1936 recvdispls(procIndex) = recvdispls(procIndex-1) + numHeaderMax
1937 end do
1938
1939 ! tasks send only from those with data
1940 if (kIndex <= stateVector_VarsLevs%mykEnd) then
1941 do procIndex = 1, mmpi_nprocs
1942 sendsizes(procIndex) = allNumHeader(procIndex)
1943 end do
1944 else
1945 sendsizes(:) = 0
1946 end if
1947 senddispls(1) = 0
1948 do procIndex = 2, mmpi_nprocs
1949 senddispls(procIndex) = senddispls(procIndex-1) + numHeaderMax
1950 end do
1951
1952 if(mmpi_nprocs > 1) then
1953 call mpi_alltoallv(cols_send, sendsizes, senddispls, mmpi_datyp_real8, &
1954 cols_recv, recvsizes, recvdispls, mmpi_datyp_real8, &
1955 mmpi_comm_grid, ierr)
1956 else
1957 cols_recv(:,1) = cols_send(:,1)
1958 end if
1959 call utl_tmg_stop(36)
1960
1961 ! reorganize ensemble of distributed columns
1962 !$OMP PARALLEL DO PRIVATE (procIndex, kIndex2, headerIndex, headerIndex2)
1963 proc_loop: do procIndex = 1, mmpi_nprocs
1964 kIndex2 = stateVector_VarsLevs%allkBeg(procIndex) + kCount - 1
1965 if ( kIndex2 > stateVector_VarsLevs%allkEnd(procIndex) ) cycle proc_loop
1966 do headerIndex = 1, numHeader
1967 headerIndex2 = headerIndex + headerIndexBeg - 1
1968 allCols_ptr(kIndex2,headerIndex2) = cols_recv(headerIndex,procIndex)
1969 end do
1970 end do proc_loop
1971 !$OMP END PARALLEL DO
1972
1973 end do k_loop
1974
1975 ! impose a lower limit on HU
1976 if(col_varExist(column,'HU')) then
1977 do headerIndex = headerIndexBeg, headerIndexEnd
1978 column_ptr => col_getColumn(column,headerIndex,'HU')
1979 column_ptr(:) = max(column_ptr(:),col_rhumin)
1980 end do
1981 end if
1982
1983 ! impose a lower/upper limits on ALL cloud variables
1984 do varNameIndex = 1, vnl_numvarmaxCloud
1985 if(col_varExist(column, vnl_varNameListCloud(varNameIndex))) then
1986 do headerIndex = headerIndexBeg, headerIndexEnd
1987 column_ptr => col_getColumn(column,headerIndex, vnl_varNameListCloud(varNameIndex))
1988 column_ptr(:) = max(column_ptr(:), qlim_getMinValueCloud(vnl_varNameListCloud(varNameIndex)))
1989 column_ptr(:) = min(column_ptr(:), qlim_getMaxValueCloud(vnl_varNameListCloud(varNameIndex)))
1990 end do
1991 end if
1992 end do
1993
1994 ! Interpolate surface height separately, only exists on mpi task 0
1995 HeightSfcPresent: if ( stateVector_VarsLevs%HeightSfcPresent ) then
1996
1997 if ( mmpi_myid == 0 ) then
1998 varName = 'GZ'
1999 kIndexHeightSfc = 0
2000 step_loop_height: do stepIndex = 1, numStep
2001
2002 if ( maxval(interpInfo_nl%allNumHeaderUsed(stepIndex,:)) == 0 ) cycle step_loop_height
2003
2004 ! interpolate to the columns destined for all procs for all steps and one lev/var
2005 !$OMP PARALLEL DO PRIVATE (procIndex, yourNumHeader, ptr2d_r8)
2006 do procIndex = 1, mmpi_nprocs
2007 yourNumHeader = interpInfo_nl%allNumHeaderUsed(stepIndex,procIndex)
2008 if ( yourNumHeader > 0 ) then
2009 ptr2d_r8 => gsv_getHeightSfc(stateVector_VarsLevs)
2010 call myezsint_r8_nl( cols_hint(1:yourNumHeader,stepIndex,procIndex), &
2011 ptr2d_r8(:,:), interpInfo_nl, kIndexHeightSfc, stepIndex, procIndex )
2012 end if
2013 end do
2014 !$OMP END PARALLEL DO
2015
2016 end do step_loop_height
2017
2018 ! interpolate in time to the columns destined for all procs and one level/variable
2019 do procIndex = 1, mmpi_nprocs
2020 cols_send(:,procIndex) = 0.0d0
2021 do stepIndex = 1, numStep
2022 !$OMP PARALLEL DO PRIVATE (headerIndex, headerIndex2, headerUsedIndex)
2023 do headerUsedIndex = 1, interpInfo_nl%allNumHeaderUsed(stepIndex, procIndex)
2024 headerIndex = interpInfo_nl%stepProcData(procIndex, stepIndex)%allHeaderIndex(headerUsedIndex)
2025 ! just copy, since surface height same for all time steps
2026 headerIndex2 = headerIndex - allHeaderIndexBeg(procIndex) + 1
2027 cols_send(headerIndex2,procIndex) = cols_hint(headerUsedIndex,stepIndex,procIndex)
2028 end do
2029 !$OMP END PARALLEL DO
2030 end do
2031 end do
2032
2033 end if
2034
2035 ! mpi communication: scatter data from task 0
2036 nsize = numHeader
2037 if(mmpi_nprocs > 1) then
2038 do procIndex = 1, mmpi_nprocs
2039 displs(procIndex) = (procIndex - 1) * numHeaderMax
2040 nsizes(procIndex) = allNumHeader(procIndex)
2041 end do
2042 call rpn_comm_scatterv(cols_send, nsizes, displs, 'MPI_REAL8', &
2043 cols_recv, nsize, 'MPI_REAL8', &
2044 0, 'GRID', ierr)
2045 else
2046 cols_recv(:,1) = cols_send(:,1)
2047 end if
2048
2049 do headerIndex = headerIndexBeg, headerIndexEnd
2050 headerIndex2 = headerIndex - headerIndexBeg + 1
2051 call col_setHeightSfc(column, headerIndex, cols_recv(headerIndex2,1))
2052 end do
2053
2054 end if HeightSfcPresent
2055
2056 deallocate(cols_hint)
2057 deallocate(cols_send)
2058 deallocate(cols_recv)
2059 deallocate(cols_send_1proc)
2060
2061 if ( dealloc ) call s2c_deallocInterpInfo( inputStateVectorType='nl' )
2062
2063 end do OBSBATCH
2064
2065 if ( dealloc) call gsv_deallocate( statevector_VarsLevs )
2066
2067 if (slantPath_TO_nl) call pressureProfileMonotonicityCheck(obsSpaceData, column)
2068
2069 firstCall = .false.
2070
2071 if ( .not. beSilent ) then
2072 write(*,*) 's2c_nl: FINISHED'
2073 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2074 end if
2075
2076 call utl_tmg_stop(34)
2077
2078 call utl_tmg_stop(30)
2079
2080 end subroutine s2c_nl
2081
2082 ! -------------------------------------------------
2083 ! myezsint_nl: Scalar field horizontal interpolation
2084 ! -------------------------------------------------
2085 subroutine myezsint_nl( column_out, field_in, interpInfo, kIndex, stepIndex, procIndex )
2086 !
2087 !:Purpose: Scalar horizontal interpolation, replaces the
2088 ! ezsint routine from rmnlib.
2089 !
2090 implicit none
2091
2092 ! Arguments:
2093 real(8) , intent(out) :: column_out(:)
2094 real(4) , intent(in) :: field_in(:,:)
2095 type(struct_interpInfo), intent(in) :: interpInfo
2096 integer , intent(in) :: stepIndex
2097 integer , intent(in) :: procIndex
2098 integer , intent(in) :: kIndex
2099
2100 ! Locals:
2101 integer :: lonIndex, latIndex, gridptIndex, headerIndex, subGridIndex, numColumn
2102 real(8) :: interpValue, weight
2103
2104 numColumn = size( column_out )
2105
2106 do headerIndex = 1, numColumn
2107
2108 ! Interpolate the model state to the obs point
2109 interpValue = 0.0d0
2110
2111 do subGridIndex = 1, interpInfo%hco%numSubGrid
2112
2113 do gridptIndex = &
2114 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex), &
2115 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex, headerIndex, kIndex)
2116
2117 lonIndex = interpInfo%lonIndexDepot(gridptIndex)
2118 latIndex = interpInfo%latIndexDepot(gridptIndex)
2119 weight = interpInfo%interpWeightDepot(gridptIndex)
2120
2121 interpValue = interpValue + weight * real(field_in(lonIndex, latIndex),8)
2122
2123 end do
2124
2125 end do
2126 column_out(headerIndex) = interpValue
2127
2128 end do
2129
2130 end subroutine myezsint_nl
2131
2132 ! -------------------------------------------------
2133 ! myezsint_r8_nl: Scalar field horizontal interpolation
2134 ! -------------------------------------------------
2135 subroutine myezsint_r8_nl( column_out, field_in, interpInfo, kIndex, stepIndex, procIndex )
2136 !
2137 !:Purpose: Scalar horizontal interpolation, replaces the
2138 ! ezsint routine from rmnlib.
2139 !
2140 implicit none
2141
2142 ! Arguments:
2143 real(8) , intent(out) :: column_out(:)
2144 real(8) , intent(in) :: field_in(:,:)
2145 type(struct_interpInfo), intent(in) :: interpInfo
2146 integer , intent(in) :: stepIndex
2147 integer , intent(in) :: procIndex
2148 integer , intent(in) :: kIndex
2149
2150 ! Locals:
2151 integer :: lonIndex, latIndex, gridptIndex, headerIndex, subGridIndex, numColumn
2152 real(8) :: interpValue, weight
2153
2154 numColumn = size( column_out )
2155
2156 do headerIndex = 1, numColumn
2157
2158 ! Interpolate the model state to the obs point
2159 interpValue = 0.0d0
2160
2161 do subGridIndex = 1, interpInfo%hco%numSubGrid
2162
2163 do gridptIndex = &
2164 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex), &
2165 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex, headerIndex, kIndex)
2166
2167 lonIndex = interpInfo%lonIndexDepot(gridptIndex)
2168 latIndex = interpInfo%latIndexDepot(gridptIndex)
2169 weight = interpInfo%interpWeightDepot(gridptIndex)
2170
2171 interpValue = interpValue + weight * field_in(lonIndex, latIndex)
2172
2173 end do
2174
2175 end do
2176 column_out(headerIndex) = interpValue
2177
2178 end do
2179
2180 end subroutine myezsint_r8_nl
2181
2182 ! -------------------------------------------------
2183 ! myezsint_tl: Scalar field horizontal interpolation
2184 ! -------------------------------------------------
2185 subroutine myezsint_tl( column_out, field_in, interpInfo, kIndex, stepIndex, procIndex )
2186 !
2187 !:Purpose: Scalar horizontal interpolation, replaces the
2188 ! ezsint routine from rmnlib.
2189 !
2190 implicit none
2191
2192 ! Arguments:
2193 real(8) , intent(out) :: column_out(:)
2194 real(pre_incrReal) , intent(in) :: field_in(:,:)
2195 type(struct_interpInfo), intent(in) :: interpInfo
2196 integer , intent(in) :: stepIndex
2197 integer , intent(in) :: procIndex
2198 integer , intent(in) :: kIndex
2199
2200 ! Locals:
2201 integer :: lonIndex, latIndex, gridptIndex, headerIndex, subGridIndex, numColumn
2202 real(8) :: interpValue, weight
2203
2204 numColumn = size( column_out )
2205
2206 do headerIndex = 1, numColumn
2207
2208 ! Interpolate the model state to the obs point
2209 interpValue = 0.0d0
2210
2211 do subGridIndex = 1, interpInfo%hco%numSubGrid
2212
2213 do gridptIndex = &
2214 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex), &
2215 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex, headerIndex, kIndex)
2216
2217 lonIndex = interpInfo%lonIndexDepot(gridptIndex)
2218 latIndex = interpInfo%latIndexDepot(gridptIndex)
2219 weight = interpInfo%interpWeightDepot(gridptIndex)
2220
2221 interpValue = interpValue + weight * field_in(lonIndex, latIndex)
2222
2223 end do
2224
2225 end do
2226 column_out(headerIndex) = interpValue
2227
2228 end do
2229
2230 end subroutine myezsint_tl
2231
2232 ! -------------------------------------------------------------
2233 ! myezsint_ad: Adjoint of scalar field horizontal interpolation
2234 ! -------------------------------------------------------------
2235 subroutine myezsint_ad( column_in, field_out, interpInfo, kIndex, stepIndex, procIndex )
2236 !
2237 !:Purpose: Adjoint of the scalar horizontal interpolation.
2238 !
2239 implicit none
2240
2241 ! Arguments:
2242 real(8) , intent(in) :: column_in(:)
2243 real(pre_incrReal) , intent(inout) :: field_out(:,:)
2244 type(struct_interpInfo), intent(in) :: interpInfo
2245 integer , intent(in) :: stepIndex
2246 integer , intent(in) :: procIndex
2247 integer , intent(in) :: kIndex
2248
2249 ! Locals:
2250 integer :: lonIndex, latIndex, gridptIndex, headerIndex, subGridIndex, numColumn
2251 real(8) :: weight
2252
2253 numColumn = size( column_in )
2254
2255 do headerIndex = 1, numColumn
2256
2257 ! Interpolate the model state to the obs point
2258
2259 do subGridIndex = 1, interpInfo%hco%numSubGrid
2260
2261 do gridptIndex = &
2262 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex), &
2263 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex, headerIndex, kIndex)
2264
2265 lonIndex = interpInfo%lonIndexDepot(gridptIndex)
2266 latIndex = interpInfo%latIndexDepot(gridptIndex)
2267 weight = interpInfo%interpWeightDepot(gridptIndex)
2268
2269 field_out(lonIndex, latIndex) = field_out(lonIndex, latIndex) + &
2270 weight * column_in(headerIndex)
2271
2272 end do
2273
2274 end do
2275
2276 end do
2277
2278 end subroutine myezsint_ad
2279
2280 ! -------------------------------------------------------------
2281 ! myezuvint_nl: Vector field horizontal interpolation
2282 ! -------------------------------------------------------------
2283 subroutine myezuvint_nl( column_out, varName, fieldUU_in, fieldVV_in, &
2284 interpInfo, kIndex, stepIndex, procIndex )
2285 !
2286 !:Purpose: Vector horizontal interpolation, replaces the
2287 ! ezuvint routine from rmnlib.
2288 !
2289 implicit none
2290
2291 ! Arguments:
2292 real(8) , intent(out) :: column_out(:)
2293 character(len=*) , intent(in) :: varName
2294 real(4) , intent(in) :: fieldUU_in(:,:)
2295 real(4) , intent(in) :: fieldVV_in(:,:)
2296 type(struct_interpInfo), intent(in) :: interpInfo
2297 integer , intent(in) :: stepIndex
2298 integer , intent(in) :: procIndex
2299 integer , intent(in) :: kIndex
2300
2301 ! Locals:
2302 integer :: lonIndex, latIndex, indexBeg, indexEnd, gridptIndex, headerIndex
2303 integer :: numColumn, subGridIndex
2304 real(8) :: interpUU(interpInfo%hco%numSubGrid), interpVV(interpInfo%hco%numSubGrid)
2305 real(8) :: lat, lon, latRot, lonRot, weight
2306 logical :: doUU, doVV
2307
2308 numColumn = size( column_out )
2309
2310 doUU = (trim(varName) == 'UU' .or. interpInfo%hco%rotated)
2311 doVV = (trim(varName) == 'VV' .or. interpInfo%hco%rotated)
2312
2313 header_loop: do headerIndex = 1, numColumn
2314
2315 interpUU(:) = 0.0d0
2316 interpVV(:) = 0.0d0
2317
2318 subGrid_loop: do subGridIndex = 1, interpInfo%hco%numSubGrid
2319
2320 indexBeg = interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex)
2321 indexEnd = interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex, headerIndex, kIndex)
2322
2323 if ( indexEnd < IndexBeg ) cycle subGrid_loop
2324
2325 ! Interpolate the model UU to the obs point
2326 do gridptIndex = indexBeg, indexEnd
2327
2328 lonIndex = interpInfo%lonIndexDepot(gridptIndex)
2329 latIndex = interpInfo%latIndexDepot(gridptIndex)
2330 weight = interpInfo%interpWeightDepot(gridptIndex)
2331
2332 if ( doUU ) interpUU(subGridIndex) = interpUU(subGridIndex) + &
2333 weight * real(fieldUU_in(lonIndex, latIndex),8)
2334 if ( doVV ) interpVV(subGridIndex) = interpVV(subGridIndex) + &
2335 weight * real(fieldVV_in(lonIndex, latIndex),8)
2336
2337 end do
2338 ! now rotate the wind vector
2339 if ( interpInfo%hco%rotated ) then
2340 lat = interpInfo%stepProcData(procIndex, stepIndex)%allLat(headerIndex, kIndex)
2341 lon = interpInfo%stepProcData(procIndex, stepIndex)%allLon(headerIndex, kIndex)
2342 latRot = interpInfo%stepProcData(procIndex,stepIndex)%allLatRot(subGridIndex, headerIndex, kIndex)
2343 lonRot = interpInfo%stepProcData(procIndex,stepIndex)%allLonRot(subGridIndex, headerIndex, kIndex)
2344
2345 call uvr_rotateWind_nl( interpInfo%uvr, & ! IN
2346 subGridIndex, & ! IN
2347 interpUU(subGridIndex), & ! INOUT
2348 interpVV(subGridIndex), & ! INOUT
2349 lat, lon, latRot, lonRot, & ! IN
2350 'ToMetWind' ) ! IN
2351 end if
2352
2353 end do subGrid_loop
2354
2355 ! return only the desired component
2356 if ( trim(varName) == 'UU' ) then
2357 column_out(headerIndex) = sum(interpUU(:))
2358 else
2359 column_out(headerIndex) = sum(interpVV(:))
2360 end if
2361
2362 end do header_loop
2363
2364 end subroutine myezuvint_nl
2365
2366 ! -------------------------------------------------------------
2367 ! myezuvint_tl: Vector field horizontal interpolation
2368 ! -------------------------------------------------------------
2369 subroutine myezuvint_tl( column_out, varName, fieldUU_in, fieldVV_in, &
2370 interpInfo, kIndex, stepIndex, procIndex )
2371 !:Purpose: Vector horizontal interpolation, replaces the
2372 ! ezuvint routine from rmnlib.
2373 !
2374 implicit none
2375
2376 ! Arguments:
2377 real(8) , intent(out) :: column_out(:)
2378 character(len=*) , intent(in) :: varName
2379 real(pre_incrReal) , intent(in) :: fieldUU_in(:,:)
2380 real(pre_incrReal) , intent(in) :: fieldVV_in(:,:)
2381 type(struct_interpInfo), intent(in) :: interpInfo
2382 integer , intent(in) :: stepIndex
2383 integer , intent(in) :: procIndex
2384 integer , intent(in) :: kIndex
2385
2386 ! Locals:
2387 integer :: lonIndex, latIndex, indexBeg, indexEnd, gridptIndex, headerIndex
2388 integer :: numColumn, subGridIndex
2389 real(8) :: interpUU(interpInfo%hco%numSubGrid), interpVV(interpInfo%hco%numSubGrid)
2390 real(8) :: lat, lon, latRot, lonRot, weight
2391 logical :: doUU, doVV
2392
2393 numColumn = size( column_out )
2394
2395 doUU = (trim(varName) == 'UU' .or. interpInfo%hco%rotated)
2396 doVV = (trim(varName) == 'VV' .or. interpInfo%hco%rotated)
2397
2398 header_loop: do headerIndex = 1, numColumn
2399
2400 interpUU(:) = 0.0d0
2401 interpVV(:) = 0.0d0
2402
2403 subGrid_loop: do subGridIndex = 1, interpInfo%hco%numSubGrid
2404
2405 indexBeg = interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex)
2406 indexEnd = interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex, headerIndex, kIndex)
2407
2408 if ( indexEnd < IndexBeg ) cycle subGrid_loop
2409
2410 ! Interpolate the model UU to the obs point
2411 do gridptIndex = indexBeg, indexEnd
2412
2413 lonIndex = interpInfo%lonIndexDepot(gridptIndex)
2414 latIndex = interpInfo%latIndexDepot(gridptIndex)
2415 weight = interpInfo%interpWeightDepot(gridptIndex)
2416
2417 if ( doUU ) interpUU(subGridIndex) = interpUU(subGridIndex) + &
2418 weight * fieldUU_in(lonIndex, latIndex)
2419 if ( doVV ) interpVV(subGridIndex) = interpVV(subGridIndex) + &
2420 weight * fieldVV_in(lonIndex, latIndex)
2421
2422 end do
2423 ! now rotate the wind vector
2424 if ( interpInfo%hco%rotated ) then
2425 lat = interpInfo%stepProcData(procIndex, stepIndex)%allLat(headerIndex, kIndex)
2426 lon = interpInfo%stepProcData(procIndex, stepIndex)%allLon(headerIndex, kIndex)
2427 latRot = interpInfo%stepProcData(procIndex,stepIndex)%allLatRot(subGridIndex, headerIndex, kIndex)
2428 lonRot = interpInfo%stepProcData(procIndex,stepIndex)%allLonRot(subGridIndex, headerIndex, kIndex)
2429
2430 call uvr_rotateWind_tl( interpInfo%uvr, & ! IN
2431 subGridIndex, & ! IN
2432 interpUU(subGridIndex), & ! INOUT
2433 interpVV(subGridIndex), & ! INOUT
2434 lat, lon, latRot, lonRot, & ! IN
2435 'ToMetWind' ) ! IN
2436 end if
2437
2438 end do subGrid_loop
2439
2440 ! return only the desired component
2441 if ( trim(varName) == 'UU' ) then
2442 column_out(headerIndex) = sum(interpUU(:))
2443 else
2444 column_out(headerIndex) = sum(interpVV(:))
2445 end if
2446
2447 end do header_loop
2448
2449 end subroutine myezuvint_tl
2450
2451 ! -------------------------------------------------------------
2452 ! myezuvint_ad: Adjoint of vector field horizontal interpolation
2453 ! -------------------------------------------------------------
2454 subroutine myezuvint_ad( column_in, varName, fieldUU_out, fieldVV_out, &
2455 interpInfo, kIndex, stepIndex, procIndex )
2456 !
2457 !:Purpose: Adjoint of the vector horizontal interpolation.
2458 !
2459 implicit none
2460
2461 ! Arguments:
2462 real(8) , intent(in) :: column_in(:)
2463 character(len=*) , intent(in) :: varName
2464 real(pre_incrReal) , intent(inout) :: fieldUU_out(:,:)
2465 real(pre_incrReal) , intent(inout) :: fieldVV_out(:,:)
2466 type(struct_interpInfo), intent(in) :: interpInfo
2467 integer , intent(in) :: stepIndex
2468 integer , intent(in) :: procIndex
2469 integer , intent(in) :: kIndex
2470
2471 ! Locals:
2472 integer :: lonIndex, latIndex, indexBeg, indexEnd, gridptIndex, headerIndex
2473 integer :: numColumn, subGridIndex
2474 real(8) :: interpUU(interpInfo%hco%numSubGrid), interpVV(interpInfo%hco%numSubGrid)
2475 real(8) :: lat, lon, latRot, lonRot, weight
2476 logical :: doUU, doVV
2477
2478 numColumn = size( column_in )
2479
2480 doUU = (trim(varName) == 'UU' .or. interpInfo%hco%rotated)
2481 doVV = (trim(varName) == 'VV' .or. interpInfo%hco%rotated)
2482
2483 header_loop: do headerIndex = 1, numColumn
2484
2485 if ( trim(varName) == 'UU' ) then
2486 interpUU(:) = column_in(headerIndex)
2487 interpVV(:) = 0.0d0
2488 else
2489 interpUU(:) = 0.0d0
2490 interpVV(:) = column_in(headerIndex)
2491 end if
2492
2493 subGrid_loop: do subGridIndex = 1, interpInfo%hco%numSubGrid
2494
2495 indexBeg = interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex)
2496 indexEnd = interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex, headerIndex, kIndex)
2497
2498 if ( indexEnd < IndexBeg ) cycle subGrid_loop
2499
2500 ! now rotate the wind vector and return the desired component
2501 if ( interpInfo%hco%rotated ) then
2502 lat = interpInfo%stepProcData(procIndex, stepIndex)%allLat(headerIndex, kIndex)
2503 lon = interpInfo%stepProcData(procIndex, stepIndex)%allLon(headerIndex, kIndex)
2504 latRot = interpInfo%stepProcData(procIndex,stepIndex)%allLatRot(subGridIndex, headerIndex, kIndex)
2505 lonRot = interpInfo%stepProcData(procIndex,stepIndex)%allLonRot(subGridIndex, headerIndex, kIndex)
2506
2507 call uvr_rotateWind_ad( interpInfo%uvr, & ! IN
2508 subGridIndex, & ! IN
2509 interpUU(subGridIndex), & ! INOUT
2510 interpVV(subGridIndex), & ! INOUT
2511 lat, lon, latRot, lonRot, & ! IN
2512 'ToMetWind' ) ! IN
2513 end if
2514
2515 ! Interpolate the model VV to the obs point
2516 do gridptIndex = indexBeg, indexEnd
2517
2518 lonIndex = interpInfo%lonIndexDepot(gridptIndex)
2519 latIndex = interpInfo%latIndexDepot(gridptIndex)
2520 weight = interpInfo%interpWeightDepot(gridptIndex)
2521
2522 if ( doUU ) fieldUU_out(lonIndex, latIndex) = &
2523 fieldUU_out(lonIndex, latIndex) + weight * interpUU(subGridIndex)
2524 if ( doVV ) fieldVV_out(lonIndex, latIndex) = &
2525 fieldVV_out(lonIndex, latIndex) + weight * interpVV(subGridIndex)
2526
2527 end do
2528
2529 end do subGrid_loop
2530
2531 end do header_loop
2532
2533 end subroutine myezuvint_ad
2534
2535 !---------------------------------------------------------
2536 ! s2c_bgcheck_bilin
2537 !---------------------------------------------------------
2538 subroutine s2c_bgcheck_bilin(column,statevector,obsSpaceData)
2539 !
2540 !:Purpose: Special version of s2c_tl used for background check. This should
2541 ! be replaced by direct call to s2c_tl. It is not general enough to
2542 ! be used for new analysis variables.
2543 !
2544 implicit none
2545
2546 ! Arguments:
2547 type(struct_columnData), intent(inout) :: column
2548 type(struct_gsv), intent(in) :: statevector
2549 type(struct_obs), intent(in) :: obsSpaceData
2550
2551 ! Locals:
2552 integer :: jk, jk2, jgl, headerIndex
2553 integer :: lonIndex, ila, ierr, subGridIndex
2554 integer :: extraLongitude
2555 real(8) :: lat, lon
2556 real(4) :: lat_r4, lon_r4, lat_deg_r4, lon_deg_r4, xpos_r4, ypos_r4, xpos2_r4, ypos2_r4
2557 real(8) :: dldy, dlw1, dlw2, dlw3, dlw4, dldx, ypos, xpos
2558 real(8), allocatable ::zgd(:,:,:)
2559 real(8), pointer :: field_ptr(:,:,:)
2560 real(8), pointer :: varColumn(:)
2561 integer :: varIndex
2562 character(len=4) :: varName
2563
2564 call utl_tmg_start(30,'--StateToColumn')
2565
2566 ! Note: We assume here the all the obs between the poles and the last grid points
2567 ! (i.e. outside the grid) have been moved within the grid by suprep
2568
2569 if (statevector%hco%global) then
2570 extraLongitude = 1
2571 else
2572 extraLongitude = 0
2573 end if
2574
2575 allocate(zgd(statevector%ni+extraLongitude,statevector%nj,statevector%nk))
2576
2577 zgd(:,:,:)=0.0d0
2578 call gsv_getField(statevector,field_ptr)
2579 zgd(1:statevector%ni,1:statevector%nj,1:statevector%nk)= &
2580 field_ptr(1:statevector%ni,1:statevector%nj,1:statevector%nk)
2581
2582 !
2583 !- 1. Expand field by repeating meridian 1 into into meridian ni+1
2584 !
2585 if (extraLongitude == 1) then
2586 do jk = 1, statevector%nk
2587 do jgl = 1, statevector%nj
2588 zgd(statevector%ni+1,jgl,jk) = zgd( 1,jgl,jk)
2589 end do
2590 end do
2591 end if
2592
2593 !
2594 !- 2. Loop over all the headers
2595 !
2596 do headerIndex = 1, col_getNumCol(column)
2597
2598 !- 2.1 Find the obs positin within the analysis grid
2599 lat = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
2600 lon = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
2601 lat_r4 = real(lat,4)
2602 lon_r4 = real(lon,4)
2603 if (lon_r4.lt.0.0 ) lon_r4 = lon_r4 + 2.0*MPC_PI_R4
2604 if (lon_r4.ge.2.*MPC_PI_R4) lon_r4 = lon_r4 - 2.0*MPC_PI_R4
2605 lat_deg_r4 = lat_r4 * MPC_DEGREES_PER_RADIAN_R4 ! Radian To Degree
2606 lon_deg_r4 = lon_r4 * MPC_DEGREES_PER_RADIAN_R4
2607 ierr = gpos_getPositionXY( stateVector % hco % EZscintID, &
2608 xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
2609 lat_deg_r4, lon_deg_r4, subGridIndex )
2610 xpos = real(xpos_r4,8)
2611 ypos = real(ypos_r4,8)
2612
2613 !- Make sure we are within bounds
2614 if ( ypos < 1.d0 .or. &
2615 ypos > real(statevector%nj , 8) .or. &
2616 xpos < 1.d0 .or. &
2617 xpos > real(statevector%ni + extraLongitude, 8) ) then
2618 write(*,*) 's2c_bgcheck_bilin: Obs outside local domain for headerIndex = ', &
2619 headerIndex
2620 write(*,*) ' obs lat, lon position = ', &
2621 Lat*MPC_DEGREES_PER_RADIAN_R8, Lon*MPC_DEGREES_PER_RADIAN_R8
2622 write(*,*) ' obs x, y position = ', &
2623 xpos, ypos
2624 write(*,*) ' domain x_end, y_end bounds = ', &
2625 statevector%ni + extraLongitude, statevector%nj
2626 call utl_abort('s2c_bgcheck_bilin')
2627 end if
2628
2629 !- 2.2 Find the lower-left grid point next to the observation
2630 if ( xpos == real(statevector%ni + extraLongitude,8) ) then
2631 lonIndex = floor(xpos) - 1
2632 else
2633 lonIndex = floor(xpos)
2634 end if
2635
2636 if ( ypos == real(statevector%nj,8) ) then
2637 ILA = floor(ypos) - 1
2638 else
2639 ILA = floor(ypos)
2640 end if
2641
2642 !- 2.3 Compute the 4 weights of the bilinear interpolation
2643 dldx = xpos - real(lonIndex,8)
2644 dldy = ypos - real(ILA,8)
2645
2646 dlw1 = (1.d0-dldx) * (1.d0-dldy)
2647 dlw2 = dldx * (1.d0-dldy)
2648 dlw3 = (1.d0-dldx) * dldy
2649 dlw4 = dldx * dldy
2650
2651 !- 2.4 Interpolate the model state to the obs point
2652
2653 do varIndex = 1, vnl_numvarmax
2654 if (.not. col_varExist(column,trim(vnl_varNameList(varIndex)))) cycle
2655 varName=trim(vnl_varNameList(varIndex))
2656 varColumn => col_getColumn(column,headerIndex,varName)
2657
2658 if(gsv_varExist(statevector,varName)) then
2659 do jk = 1, gsv_getNumLevFromVarName(statevector,varName)
2660 jk2=jk+gsv_getOffsetFromVarName(statevector,varName)
2661 varColumn(jk) = dlw1*zgd(lonIndex ,ila,jk2) &
2662 + dlw2*zgd(lonIndex+1,ila,jk2) &
2663 + dlw3*zgd(lonIndex ,ila+1,jk2) &
2664 + dlw4*zgd(lonIndex+1,ila+1,jk2)
2665 end do
2666 end if
2667
2668 nullify(varColumn)
2669 end do
2670
2671 end do
2672
2673 deallocate(zgd)
2674
2675 call utl_tmg_stop(30)
2676
2677 end subroutine s2c_bgcheck_bilin
2678
2679 !--------------------------------------------------------------------------
2680 ! s2c_setupHorizInterp
2681 !--------------------------------------------------------------------------
2682 subroutine s2c_setupHorizInterp(footprintRadius_r4, interpInfo, &
2683 stateVector, headerIndex, kIndex, stepIndex, &
2684 procIndex, numGridpt)
2685 !
2686 !:Purpose: To identify the appropriate horizontal interpolation scheme based
2687 ! on footprint radius value. Then to call the corresponding
2688 ! subroutine to determine the grid points and their associated
2689 ! weights.
2690 !
2691 implicit none
2692
2693 ! Arguments:
2694 real(4) , intent(in) :: footprintRadius_r4 ! (metres)
2695 type(struct_interpInfo), intent(inout) :: interpInfo
2696 type(struct_gsv) , intent(in) :: stateVector
2697 integer , intent(in) :: headerIndex, kIndex, stepIndex
2698 integer , intent(in) :: procIndex
2699 integer , intent(out) :: numGridpt(interpInfo%hco%numSubGrid)
2700
2701 if ( footprintRadius_r4 > 0.0 ) then
2702
2703 call s2c_setupFootprintInterp(footprintRadius_r4, interpInfo, stateVector, headerIndex, kIndex, stepIndex, procIndex, numGridpt)
2704
2705 else if ( footprintRadius_r4 == bilinearFootprint ) then
2706
2707 call s2c_setupBilinearInterp(interpInfo, stateVector, headerIndex, kIndex, stepIndex, procIndex, numGridpt)
2708
2709 else if ( footprintRadius_r4 == lakeFootprint ) then
2710
2711 call s2c_setupLakeInterp(interpInfo, stateVector, headerIndex, kIndex, stepIndex, procIndex, numGridpt)
2712
2713 else if ( footprintRadius_r4 == nearestNeighbourFootprint ) then
2714
2715 call s2c_setupNearestNeighbor(interpInfo, stateVector, headerIndex, kIndex, stepIndex, procIndex, numGridpt)
2716
2717 else
2718
2719 write(*,*) 'footprint radius = ',footprintRadius_r4
2720 call utl_abort('s2c_setupHorizInterp: footprint radius not permitted')
2721
2722 end if
2723
2724 end subroutine s2c_setupHorizInterp
2725
2726 !------------------------------------------------------------------
2727 ! s2c_getFootprintRadius
2728 !------------------------------------------------------------------
2729 function s2c_getFootprintRadius( obsSpaceData, stateVector, headerIndex ) result(fpr)
2730 !
2731 !:Purpose: To determine the footprint radius (metres) of the observation.
2732 ! In the case of bilinear horizontal interpolation,
2733 ! the returned footprint is zero (default).
2734 !
2735 implicit none
2736
2737 ! Arguments:
2738 type(struct_obs), intent(in) :: obsSpaceData
2739 type(struct_gsv), intent(in) :: stateVector
2740 integer , intent(in) :: headerIndex
2741 ! Result:
2742 real(4) :: fpr
2743
2744 ! Locals:
2745 character(len=2) :: obsFamily
2746 character(len=12) :: cstnid
2747 integer :: codeType
2748
2749 fpr = bilinearFootprint
2750
2751 obsFamily = obs_getFamily ( obsSpaceData, headerIndex_opt=headerIndex )
2752 if ( obsFamily == 'GL' ) then
2753
2754 cstnid = obs_elem_c ( obsSpaceData, 'STID' , headerIndex )
2755 codeType = obs_headElem_i( obsSpaceData, OBS_ITY, headerIndex )
2756
2757 if (index(cstnid,'DMSP') == 1) then
2758
2759 select case(cstnid)
2760 case('DMSP15')
2761 fpr = 27.5e3
2762 case('DMSP16','DMSP17','DMSP18')
2763 fpr = 29.0e3
2764 case DEFAULT
2765 call utl_abort('s2c_getFootprintRadius: UNKNOWN station id: '//cstnid)
2766 end select
2767
2768 else if (cstnid == 'GCOM-W1') then
2769
2770 fpr = 11.0e3
2771
2772 else if (cstnid(1:6) == 'METOP-') then
2773
2774 fpr = 25.0e3
2775
2776 else if (cstnid == 'noaa-19') then
2777
2778 fpr = 2.75e3
2779
2780 else if (cstnid == 'CIS_DAILY') then
2781
2782 fpr = bilinearFootprint
2783
2784 else if (cstnid == 'RS1_IMG') then
2785
2786 fpr = bilinearFootprint
2787
2788 else if (codtyp_get_name(codeType) == 'iceclake') then
2789
2790 fpr = lakeFootprint
2791
2792 else if (cstnid == 'CIS_REGIONAL') then
2793
2794 fpr = bilinearFootprint
2795
2796 else if (cstnid(1:3) == 'RCM') then
2797
2798 fpr = 0.8e3
2799
2800 else
2801
2802 call utl_abort('s2c_getFootprintRadius: UNKNOWN station id: '//cstnid)
2803
2804 end if
2805
2806 else if (obsFamily == 'HY') then
2807
2808 fpr = nearestNeighbourFootprint
2809
2810 else if (obsFamily == 'TO' .and. useFootprintForTovs ) then
2811
2812 fpr = getTovsFootprintRadius(obsSpaceData, headerIndex, beSilent_opt=.true.)
2813
2814 ! As safety margin, add 10% to maxGridSpacing before comparing to the footprint radius.
2815 if ( fpr < 1.1 * real(stateVector%hco%maxGridSpacing,4) ) fpr = bilinearFootprint
2816
2817 else
2818
2819 fpr = bilinearFootprint
2820
2821 end if
2822
2823 end function s2c_getFootprintRadius
2824
2825 !--------------------------------------------------------------------------
2826 ! s2c_rejectZeroWeightObs
2827 !--------------------------------------------------------------------------
2828 subroutine s2c_rejectZeroWeightObs(interpInfo, obsSpaceData, mykBeg, mykEnd)
2829 !
2830 !:Purpose: To flag an observation in obsSpaceData as being rejected if
2831 ! it has zero interpolation weight (usually because an ocean
2832 ! obs is touching land) on any mpi task.
2833 !
2834 implicit none
2835
2836 ! Arguments:
2837 type(struct_interpInfo), intent(inout) :: interpInfo
2838 type(struct_obs) , intent(inout) :: obsSpaceData
2839 integer , intent(in) :: mykBeg
2840 integer , intent(in) :: mykEnd
2841
2842 ! Locals:
2843 integer :: numStep, procIndex, stepIndex, headerUsedIndex, headerIndex, kIndex
2844 integer :: numHeader, numHeaderMax, bodyIndexBeg, bodyIndexEnd, bodyIndex
2845 integer :: subGridIndex, gridptIndex, ierr, nsize
2846 integer, save :: numWrites = 0
2847 logical, allocatable :: allRejectObs(:,:), allRejectObsMpiGlobal(:,:)
2848
2849 write(*,*) 's2c_rejectZeroWeightObs: Starting'
2850
2851 numHeader = obs_numheader(obsSpaceData)
2852 call rpn_comm_allreduce(numHeader, numHeaderMax, 1, &
2853 'MPI_INTEGER', 'MPI_MAX', 'GRID', ierr)
2854
2855 allocate(allRejectObs(numHeaderMax,mmpi_nprocs))
2856 allocate(allRejectObsMpiGlobal(numHeaderMax,mmpi_nprocs))
2857 allRejectObs(:,:) = .false.
2858 allRejectObsMpiGlobal(:,:) = .false.
2859
2860 numStep = size(interpInfo%stepProcData(1,:))
2861 do procIndex = 1, mmpi_nprocs
2862 do stepIndex = 1, numStep
2863 do headerUsedIndex = 1, interpInfo%allNumHeaderUsed(stepIndex,procIndex)
2864 headerIndex = interpInfo%stepProcData(procIndex,stepIndex)%allHeaderIndex(headerUsedIndex)
2865 do kIndex = mykBeg, mykEnd
2866 if (kIndex == mykBeg) allRejectObs(headerIndex,procIndex) = .true.
2867 do subGridIndex = 1, interpInfo%hco%numSubGrid
2868 do gridptIndex = &
2869 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerUsedIndex, kIndex), &
2870 interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex, headerUsedIndex, kIndex)
2871 if (interpInfo%interpWeightDepot(gridptIndex) > 0.0d0) then
2872 allRejectObs(headerIndex,procIndex) = .false.
2873 end if
2874 end do
2875 end do
2876 end do ! kIndex
2877 end do ! headerUsedIndex
2878 end do ! stepIndex
2879 end do ! procIndex
2880
2881 ! do global communication of reject flags
2882 nsize = numHeaderMax*mmpi_nprocs
2883 call rpn_comm_allreduce(allRejectObs,allRejectObsMpiGlobal,nsize,'MPI_LOGICAL','MPI_LOR','GRID',ierr)
2884
2885 ! modify obsSpaceData based on reject flags
2886 do headerIndex = 1, obs_numHeader(obsSpaceData)
2887 if (allRejectObsMpiGlobal(headerIndex,mmpi_myid+1)) then
2888
2889 numWrites = numWrites + 1
2890 if (numWrites < maxNumWrites) then
2891 write(*,*) 's2c_rejectZeroWeightObs: Rejecting OBS with zero weight, index ', headerIndex
2892 else if (numWrites == maxNumWrites) then
2893 write(*,*) 's2c_rejectZeroWeightObs: More rejects, but reached maximum number of writes to the listing.'
2894 end if
2895
2896 bodyIndexBeg = obs_headElem_i(obsSpaceData, OBS_RLN, headerIndex)
2897 bodyIndexEnd = obs_headElem_i(obsSpaceData, OBS_NLV, headerIndex) + bodyIndexBeg -1
2898 do bodyIndex = bodyIndexBeg, bodyIndexEnd
2899 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
2900 end do
2901 call obs_headSet_i(obsSpaceData, OBS_ST1, headerIndex, &
2902 ibset( obs_headElem_i(obsSpaceData, OBS_ST1, headerIndex), 05))
2903 end if
2904 end do
2905
2906 deallocate(allRejectObs)
2907 deallocate(allRejectObsMpiGlobal)
2908
2909 write(*,*) 's2c_rejectZeroWeightObs: Finished'
2910
2911 end subroutine s2c_rejectZeroWeightObs
2912
2913 !--------------------------------------------------------------------------
2914 ! s2c_setupBilinearInterp
2915 !--------------------------------------------------------------------------
2916 subroutine s2c_setupBilinearInterp(interpInfo, stateVector, headerIndex, kIndex, &
2917 stepIndex, procIndex, numGridpt)
2918 !
2919 !:Purpose: To determine the grid points and their associated weights
2920 ! for the bilinear horizontal interpolation. If mask is present
2921 ! we currently can only handle a single 2D mask (like for sea
2922 ! ice or SST analysis). Will abort if multiple ocean levels present.
2923 !
2924 implicit none
2925
2926 ! Arguments:
2927 type(struct_interpInfo), intent(inout) :: interpInfo
2928 type(struct_gsv) , intent(in) :: stateVector
2929 integer , intent(in) :: headerIndex, kIndex, stepIndex
2930 integer , intent(in) :: procIndex
2931 integer , intent(out) :: numGridpt(interpInfo%hco%numSubGrid)
2932
2933 ! Locals:
2934 integer :: depotIndex
2935 integer :: ierr, niP1
2936 integer :: latIndex, lonIndex, latIndex2, lonIndex2, lonIndexP1
2937 integer :: subGridIndex, subGridForInterp, numSubGridsForInterp
2938 integer :: ipoint, gridptCount
2939 integer :: latIndexVec(4), lonIndexVec(4)
2940 logical :: mask(2,2)
2941 real(8) :: WeightVec(4)
2942 real(8) :: dldx, dldy
2943 real(8) :: weightsSum
2944 real(4) :: lon_deg_r4, lat_deg_r4
2945 real(4) :: xpos_r4, ypos_r4, xpos2_r4, ypos2_r4
2946 integer, parameter :: leftIndex = 1, rightIndex = 2, bottomIndex = 1, topIndex = 2
2947
2948 numGridpt(:) = 0
2949
2950 lat_deg_r4 = real(interpInfo%stepProcData(procIndex, stepIndex)%allLat(headerIndex, kIndex) * &
2951 MPC_DEGREES_PER_RADIAN_R8)
2952 lon_deg_r4 = real(interpInfo%stepProcData(procIndex, stepIndex)%allLon(headerIndex, kIndex) * &
2953 MPC_DEGREES_PER_RADIAN_R8)
2954 ierr = gpos_getPositionXY( stateVector%hco%EZscintID, &
2955 xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
2956 lat_deg_r4, lon_deg_r4, subGridIndex )
2957
2958 ! Allow for periodicity in Longitude for global Gaussian grid
2959 if ( stateVector%hco%grtyp == 'G' .or. &
2960 (stateVector%hco%grtyp == 'Z' .and. stateVector%hco%global) ) then
2961 niP1 = statevector%ni + 1
2962 else
2963 niP1 = statevector%ni
2964 end if
2965
2966 ! Find the lower-left grid point next to the observation
2967 if ( xpos_r4 >= real(niP1) ) then
2968 xpos_r4 = real(niP1)
2969 lonIndex = niP1 - 1
2970 else if ( xpos_r4 < 1.0 ) then
2971 xpos_r4 = 1.0
2972 lonIndex = 1
2973 else
2974 lonIndex = floor(xpos_r4)
2975 end if
2976 if ( xpos2_r4 >= real(niP1) ) then
2977 xpos2_r4 = real(niP1)
2978 lonIndex2 = niP1 - 1
2979 else if ( xpos2_r4 < 1.0 ) then
2980 xpos2_r4 = 1.0
2981 lonIndex2 = 1
2982 else
2983 lonIndex2 = floor(xpos2_r4)
2984 end if
2985
2986 if ( ypos_r4 >= real(statevector%nj) ) then
2987 ypos_r4 = real(statevector%nj)
2988 latIndex = statevector%nj - 1
2989 else if ( ypos_r4 < 1.0 ) then
2990 ypos_r4 = 1.0
2991 latIndex = 1
2992 else
2993 latIndex = floor(ypos_r4)
2994 end if
2995 if ( ypos2_r4 >= real(statevector%nj) ) then
2996 ypos2_r4 = real(statevector%nj)
2997 latIndex2 = statevector%nj - 1
2998 else if ( ypos2_r4 < 1.0 ) then
2999 ypos2_r4 = 1.0
3000 latIndex2 = 1
3001 else
3002 latIndex2 = floor(ypos2_r4)
3003 end if
3004
3005 if ( stateVector%hco%grtyp == 'U' ) then
3006 if ( ypos_r4 == real(stateVector%nj/2) ) then
3007 latIndex = floor(ypos_r4) - 1
3008 end if
3009 if ( ypos2_r4 == real(stateVector%nj/2) ) then
3010 latIndex2 = floor(ypos2_r4) - 1
3011 end if
3012 end if
3013
3014 ! Handle periodicity in longitude
3015 lonIndexP1 = lonIndex + 1
3016 if ( lonIndexP1 == statevector%ni + 1 ) lonIndexP1 = 1
3017
3018 ! Check if location is in between Yin and Yang (should not happen)
3019 if ( stateVector%hco%grtyp == 'U' ) then
3020 if ( ypos_r4 > real(stateVector%nj/2) .and. &
3021 ypos_r4 < real((stateVector%nj/2)+1) ) then
3022 write(*,*) 's2c_setupBilinearInterp: WARNING, obs position in between Yin and Yang!!!'
3023 write(*,*) ' xpos, ypos = ', xpos_r4, ypos_r4
3024 end if
3025 if ( ypos2_r4 > real(stateVector%nj/2) .and. &
3026 ypos2_r4 < real((stateVector%nj/2)+1) ) then
3027 write(*,*) 's2c_setupBilinearInterp: WARNING, obs position in between Yin and Yang!!!'
3028 write(*,*) ' xpos2, ypos2 = ', xpos2_r4, ypos2_r4
3029 end if
3030 end if
3031
3032 if ( subGridIndex == 3 ) then
3033 ! both subGrids involved in interpolation, so first treat subGrid 1
3034 numSubGridsForInterp = 2
3035 subGridIndex = 1
3036 else
3037 ! only 1 subGrid involved in interpolation
3038 numSubGridsForInterp = 1
3039 end if
3040
3041 if ( stateVector%oceanMask%maskPresent ) then
3042 ! abort if 3D mask is present, since we may not handle this situation correctly
3043 if ( stateVector%oceanMask%nLev > 1 ) then
3044 call utl_abort('s2c_setupBilinearInterp: 3D mask present - this case not properly handled')
3045 end if
3046 ! extract the ocean mask
3047 mask(leftIndex ,bottomIndex) = stateVector%oceanMask%mask(lonIndex ,latIndex ,1)
3048 mask(rightIndex,bottomIndex) = stateVector%oceanMask%mask(lonIndexP1,latIndex ,1)
3049 mask(leftIndex ,topIndex ) = stateVector%oceanMask%mask(lonIndex ,latIndex + 1,1)
3050 mask(rightIndex,topIndex ) = stateVector%oceanMask%mask(lonIndexP1,latIndex + 1,1)
3051 else
3052 mask(:,:) = .true.
3053 end if
3054
3055 do subGridForInterp = 1, numSubGridsForInterp
3056
3057 WeightVec(:) = 0
3058 gridptCount = 0
3059
3060 ! Compute the 4 weights of the bilinear interpolation
3061 if ( subGridForInterp == 1 ) then
3062 ! when only 1 subGrid involved, subGridIndex can be 1 or 2
3063 dldx = real(xpos_r4,8) - real(lonIndex,8)
3064 dldy = real(ypos_r4,8) - real(latIndex,8)
3065 else
3066 ! when 2 subGrids, subGridIndex is set to 1 for 1st iteration, 2 for second
3067 subGridIndex = 2
3068 lonIndex = lonIndex2
3069 latIndex = latIndex2
3070 lonIndexP1 = lonIndex2 + 1
3071 dldx = real(xpos2_r4,8) - real(lonIndex,8)
3072 dldy = real(ypos2_r4,8) - real(latIndex,8)
3073 end if
3074
3075 if ( mask(leftIndex ,bottomIndex) ) then
3076 gridptCount = gridptCount + 1
3077 latIndexVec(gridptCount) = latIndex
3078 lonIndexVec(gridptCount) = lonIndex
3079 WeightVec(gridptCount) = (1.d0-dldx) * (1.d0-dldy)
3080 end if
3081
3082 if ( mask(rightIndex,bottomIndex) ) then
3083 gridptCount = gridptCount + 1
3084 latIndexVec(gridptCount) = latIndex
3085 lonIndexVec(gridptCount) = lonIndexP1
3086 WeightVec(gridptCount) = dldx * (1.d0-dldy)
3087 end if
3088
3089 if ( mask(leftIndex ,topIndex ) ) then
3090 gridptCount = gridptCount + 1
3091 latIndexVec(gridptCount) = latIndex + 1
3092 lonIndexVec(gridptCount) = lonIndex
3093 WeightVec(gridptCount) = (1.d0-dldx) * dldy
3094 end if
3095
3096 if ( mask(rightIndex,topIndex ) ) then
3097 gridptCount = gridptCount + 1
3098 latIndexVec(gridptCount) = latIndex + 1
3099 lonIndexVec(gridptCount) = lonIndexP1
3100 WeightVec(gridptCount) = dldx * dldy
3101 end if
3102
3103 weightsSum = sum(WeightVec(1:gridptCount))
3104 if ( weightsSum > 0.d0 ) then
3105 WeightVec(1:gridptCount) = WeightVec(1:gridptCount) / weightsSum
3106 end if
3107
3108 ! divide weight by number of subGrids
3109 WeightVec(1:gridptCount) = WeightVec(1:gridptCount) / real(numSubGridsForInterp,8)
3110
3111 if ( allocated(interpInfo%interpWeightDepot) ) then
3112
3113 depotIndex = interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex)
3114
3115 do ipoint=1,gridptCount
3116
3117 interpInfo%interpWeightDepot(depotIndex) = WeightVec(ipoint)
3118 interpInfo%latIndexDepot(depotIndex) = latIndexVec(ipoint)
3119 interpInfo%lonIndexDepot(depotIndex) = lonIndexVec(ipoint)
3120 depotIndex = depotIndex + 1
3121
3122 end do
3123
3124 end if
3125
3126 numGridpt(subGridIndex) = gridptCount
3127
3128 end do ! subGrid
3129
3130 end subroutine s2c_setupBilinearInterp
3131
3132 !--------------------------------------------------------------------------
3133 ! s2c_setupFootprintInterp
3134 !--------------------------------------------------------------------------
3135 subroutine s2c_setupFootprintInterp(fpr, interpInfo, stateVector, headerIndex, &
3136 kIndex, stepIndex, procIndex, numGridpt)
3137 !
3138 !:Purpose: To determine the grid points and their associated weights
3139 ! for the footprint horizontal interpolation.
3140 !
3141 implicit none
3142
3143 ! Arguments:
3144 real(4) , intent(in) :: fpr ! footprint radius (metres)
3145 type(struct_interpInfo), intent(inout) :: interpInfo
3146 type(struct_gsv) , intent(in) :: stateVector
3147 integer , intent(in) :: headerIndex, kIndex, stepIndex
3148 integer , intent(in) :: procIndex
3149 integer , intent(out) :: numGridpt(interpInfo%hco%numSubGrid)
3150
3151 ! Locals:
3152 integer :: depotIndex
3153 integer :: ierr
3154 integer :: latIndexCentre, lonIndexCentre, latIndexCentre2, lonIndexCentre2
3155 integer :: subGridIndex, numLocalGridptsFoundSearch
3156 real(4) :: lonObs_deg_r4, latObs_deg_r4
3157 real(8) :: lonObs, latObs
3158 real(4) :: xpos_r4, ypos_r4, xpos2_r4, ypos2_r4
3159 integer :: ipoint, gridptCount
3160 integer :: lonIndex, latIndex, resultsIndex, gridIndex
3161 integer :: lonIndexVec(maxNumLocalGridptsSearch), latIndexVec(maxNumLocalGridptsSearch)
3162 type(kdtree2_result) :: searchResults(maxNumLocalGridptsSearch)
3163 real(kdkind) :: refPosition(3), maxRadiusSquared
3164 type(kdtree2), pointer :: tree
3165
3166 numGridpt(:) = 0
3167
3168 ! Determine the grid point nearest the observation.
3169
3170 latObs = interpInfo % stepProcData(procIndex, stepIndex) % allLat(headerIndex, kIndex)
3171 lonObs = interpInfo % stepProcData(procIndex, stepIndex) % allLon(headerIndex, kIndex)
3172
3173 latObs_deg_r4 = real(latObs * MPC_DEGREES_PER_RADIAN_R8)
3174 lonObs_deg_r4 = real(lonObs * MPC_DEGREES_PER_RADIAN_R8)
3175 ierr = gpos_getPositionXY( stateVector%hco%EZscintID, &
3176 xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
3177 latObs_deg_r4, lonObs_deg_r4, subGridIndex )
3178
3179 lonIndexCentre = nint(xpos_r4)
3180 latIndexCentre = nint(ypos_r4)
3181 lonIndexCentre2 = nint(xpos2_r4)
3182 latIndexCentre2 = nint(ypos2_r4)
3183
3184 if ( subGridIndex == 3 ) then
3185 write(*,*) 's2c_setupFootprintInterp: revise code'
3186 call utl_abort('s2c_setupFootprintInterp: both subGrids involved in interpolation.')
3187 end if
3188
3189 ! Return if observation is not on the grid, or masked.
3190 if ( lonIndexCentre < 1 .or. lonIndexCentre > statevector%hco%ni .or. &
3191 latIndexCentre < 1 .or. latIndexCentre > statevector%hco%nj ) return
3192
3193 if ( stateVector%oceanMask%maskPresent ) then
3194 ! abort if 3D mask is present, since we may not handle this situation correctly
3195 if ( stateVector%oceanMask%nLev > 1 ) then
3196 call utl_abort('s2c_setupFootprintInterp: 3D mask present - this case not properly handled')
3197 end if
3198
3199 if ( .not. stateVector%oceanMask%mask(lonIndexCentre,latIndexCentre,1) ) return
3200 end if
3201
3202 ! do the search
3203 maxRadiusSquared = real(fpr,8) ** 2
3204 refPosition(:) = kdtree2_3dPosition(lonObs, latObs)
3205 nullify(tree)
3206 if ( interpInfo%inputStateVectorType == 'nl' ) then
3207 if ( associated(tree_nl) ) then
3208 tree => tree_nl
3209 else
3210 call utl_abort('s2c_setupFootprintInterp: tree_nl is not allocated!')
3211 end if
3212 else if ( interpInfo%inputStateVectorType == 'tl' .or. &
3213 interpInfo%inputStateVectorType == 'ad' ) then
3214 if ( associated(tree_tlad) ) then
3215 tree => tree_tlad
3216 else
3217 call utl_abort('s2c_setupFootprintInterp: tree_tlad is not allocated!')
3218 end if
3219 end if
3220 call kdtree2_r_nearest(tp=tree, qv=refPosition, r2=maxRadiusSquared, &
3221 nfound=numLocalGridptsFoundSearch, &
3222 nalloc=maxNumLocalGridptsSearch, &
3223 results=searchResults)
3224
3225 if (numLocalGridptsFoundSearch > maxNumLocalGridptsSearch ) then
3226 call utl_abort('s2c_setupFootprintInterp: the parameter maxNumLocalGridptsSearch must be increased')
3227 else if ( numLocalGridptsFoundSearch < minNumLocalGridptsSearch .and. useFootprintForTovs ) then
3228 write(*,*) 's2c_setupFootprintInterp: Warning! For TOVS headerIndex=', headerIndex, &
3229 ' number of grid points found within footprint radius=', fpr, ' is less than ', &
3230 minNumLocalGridptsSearch
3231 end if
3232
3233 ! ensure at least the nearest neighbor is included in lonIndexVec/latIndexVec
3234 ! if footprint size is smaller than the grid spacing.
3235 gridptCount = 1
3236 lonIndexVec(gridptCount) = lonIndexCentre
3237 latIndexVec(gridptCount) = latIndexCentre
3238
3239 ! fill the rest of lonIndexVec/latIndexVec
3240 gridLoop1: do resultsIndex = 1, numLocalGridptsFoundSearch
3241 gridIndex = searchResults(resultsIndex)%idx
3242 if ( gridIndex < 1 .or. gridIndex > statevector%hco%ni * statevector%hco%nj ) then
3243 write(*,*) 's2c_setupFootprintInterp: gridIndex=', gridIndex
3244 call utl_abort('s2c_setupFootprintInterp: gridIndex out of bound.')
3245 end if
3246
3247 latIndex = (gridIndex - 1) / statevector%hco%ni + 1
3248 lonIndex = gridIndex - (latIndex - 1) * statevector%hco%ni
3249 if ( lonIndex < 1 .or. lonIndex > statevector%hco%ni .or. &
3250 latIndex < 1 .or. latIndex > statevector%hco%nj ) then
3251 write(*,*) 's2c_setupFootprintInterp: lonIndex=', lonIndex, ',latIndex=', latIndex
3252 call utl_abort('s2c_setupFootprintInterp: lonIndex/latIndex out of bound.')
3253 end if
3254
3255 if ( stateVector%oceanMask%maskPresent ) then
3256 if ( .not. stateVector%oceanMask%mask(lonIndex,latIndex,1) ) cycle gridLoop1
3257 end if
3258
3259 if ( lonIndex == lonIndexCentre .and. latIndex == latIndexCentre ) cycle gridLoop1
3260
3261 gridptCount = gridptCount + 1
3262 lonIndexVec(gridptCount) = lonIndex
3263 latIndexVec(gridptCount) = latIndex
3264 end do gridLoop1
3265
3266 if ( allocated(interpInfo%interpWeightDepot) ) then
3267
3268 depotIndex = interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex)
3269
3270 do ipoint = 1, gridptCount
3271
3272 interpInfo%interpWeightDepot(depotIndex) = 1.0d0 / real(gridptCount,8)
3273 interpInfo%latIndexDepot(depotIndex) = latIndexVec(ipoint)
3274 interpInfo%lonIndexDepot(depotIndex) = lonIndexVec(ipoint)
3275 depotIndex = depotIndex + 1
3276
3277 end do
3278
3279 end if
3280
3281 numGridpt(subGridIndex) = gridptCount
3282
3283 end subroutine s2c_setupFootprintInterp
3284
3285 !--------------------------------------------------------------------------
3286 ! s2c_setupLakeInterp
3287 !--------------------------------------------------------------------------
3288 subroutine s2c_setupLakeInterp(interpInfo, stateVector, headerIndex, kIndex, &
3289 stepIndex, procIndex, numGridpt)
3290 !
3291 !:Purpose: To determine the grid points and their associated weights
3292 ! for the lake horizontal interpolation.
3293 !
3294 implicit none
3295
3296 ! Arguments:
3297 type(struct_interpInfo), intent(inout) :: interpInfo
3298 type(struct_gsv) , intent(in) :: stateVector
3299 integer , intent(in) :: headerIndex, kIndex, stepIndex
3300 integer , intent(in) :: procIndex
3301 integer , intent(out) :: numGridpt(interpInfo%hco%numSubGrid)
3302
3303 ! Locals:
3304 integer :: depotIndex
3305 integer :: ierr
3306 integer :: latIndexCentre, lonIndexCentre, latIndexCentre2, lonIndexCentre2
3307 integer :: subGridIndex, subGridForInterp, numSubGridsForInterp
3308 real(4) :: lon_deg_r4, lat_deg_r4
3309 real(8) :: lon_rad, lat_rad
3310 real(4) :: xpos_r4, ypos_r4, xpos2_r4, ypos2_r4
3311 integer :: ipoint, gridptCount
3312 integer :: lakeCount
3313 integer :: lonIndex, latIndex, lakeIndex
3314 integer :: lonIndexVec(statevector%ni*statevector%nj), latIndexVec(statevector%ni*statevector%nj)
3315 logical :: reject, lake(statevector%ni,statevector%nj)
3316 integer :: k, l
3317
3318 if ( stateVector%hco%grtyp == 'U' ) then
3319 call utl_abort('s2c_setupLakeInterp: Yin-Yang grid not supported')
3320 end if
3321
3322 if ( .not.stateVector%oceanMask%maskPresent ) then
3323 call utl_abort('s2c_setupLakeInterp: Only compatible when mask present')
3324 end if
3325
3326 numGridpt(:) = 0
3327
3328 reject = .false.
3329
3330 numGridpt(:) = 0
3331
3332 ! Determine the grid point nearest the observation.
3333
3334 lat_rad = interpInfo%stepProcData(procIndex, stepIndex)%allLat(headerIndex, kIndex)
3335 lon_rad = interpInfo%stepProcData(procIndex, stepIndex)%allLon(headerIndex, kIndex)
3336 lat_deg_r4 = real(lat_rad * MPC_DEGREES_PER_RADIAN_R8)
3337 lon_deg_r4 = real(lon_rad * MPC_DEGREES_PER_RADIAN_R8)
3338 ierr = gpos_getPositionXY( stateVector%hco%EZscintID, &
3339 xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
3340 lat_deg_r4, lon_deg_r4, subGridIndex )
3341
3342 lonIndexCentre = nint(xpos_r4)
3343 latIndexCentre = nint(ypos_r4)
3344 lonIndexCentre2 = nint(xpos2_r4)
3345 latIndexCentre2 = nint(ypos2_r4)
3346
3347 if ( subGridIndex == 3 ) then
3348 write(*,*) 's2c_setupLakeInterp: revise code'
3349 call utl_abort('s2c_setupLakeInterp: both subGrids involved in interpolation.')
3350 numSubGridsForInterp = 2
3351 subGridIndex = 1
3352 else
3353 ! only 1 subGrid involved in interpolation
3354 numSubGridsForInterp = 1
3355 end if
3356
3357 do subGridForInterp = 1, numSubGridsForInterp
3358
3359 gridptCount = 0
3360
3361 ! It can happen that the lake location is closest to a grid point
3362 ! where MASK(I,J) = .false. while there are other grid points for the
3363 ! same lake where MASK(I,J) = .true.. Code needs modifications
3364 ! for this case.
3365
3366 ! If observation is not on the grid, don't use it.
3367 if ( lonIndexCentre < 1 .or. lonIndexCentre > statevector%ni .or. &
3368 latIndexCentre < 1 .or. latIndexCentre > statevector%nj ) reject = .true.
3369
3370 if ( .not. stateVector%oceanMask%mask(lonIndexCentre,latIndexCentre,1) ) reject = .true.
3371
3372 if ( .not. reject ) then
3373
3374 lake(:,:) = .false.
3375 lake(lonIndexCentre,latIndexCentre) = .true.
3376 gridptCount = 1
3377 lonIndexVec(gridptCount) = lonIndexCentre
3378 latIndexVec(gridptCount) = latIndexCentre
3379
3380 lakeCount = 0
3381
3382 do while(lakeCount /= gridptCount)
3383
3384 do lakeIndex = lakeCount+1, gridptCount
3385
3386 if(lakeIndex == lakeCount+1) lakeCount = gridptCount
3387
3388 k = lonIndexVec(lakeIndex)
3389 l = latIndexVec(lakeIndex)
3390
3391 do latIndex = max(1,l-1),min(l+1,statevector%nj)
3392 do lonIndex = max(1,k-1),min(k+1,statevector%ni)
3393 if(stateVector%oceanMask%mask(lonIndex,latIndex,1) .and. .not. lake(lonIndex,latIndex)) then
3394 lake(lonIndex,latIndex) = .true.
3395 gridptCount = gridptCount + 1
3396 lonIndexVec(gridptCount) = lonIndex
3397 latIndexVec(gridptCount) = latIndex
3398 end if
3399 end do
3400 end do
3401
3402 end do
3403
3404 end do
3405
3406 if ( allocated(interpInfo%interpWeightDepot) ) then
3407
3408 depotIndex = interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex)
3409
3410 do ipoint=1,gridptCount
3411
3412 interpInfo%interpWeightDepot(depotIndex) = 1.0d0 / real(gridptCount,8)
3413 interpInfo%latIndexDepot(depotIndex) = latIndexVec(ipoint)
3414 interpInfo%lonIndexDepot(depotIndex) = lonIndexVec(ipoint)
3415 depotIndex = depotIndex + 1
3416
3417 end do
3418
3419 end if
3420
3421 numGridpt(subGridIndex) = gridptCount
3422
3423 end if ! not reject
3424
3425 end do ! subGrid
3426
3427 end subroutine s2c_setupLakeInterp
3428
3429 !--------------------------------------------------------------------------
3430 ! s2c_setupNearestNeighbor
3431 !--------------------------------------------------------------------------
3432 subroutine s2c_setupNearestNeighbor(interpInfo, stateVector, headerIndex, kIndex, &
3433 stepIndex, procIndex, numGridpt)
3434 !
3435 !:Purpose: Determine the nearest grid points to the observations location
3436 !
3437 implicit none
3438
3439 ! Arguments:
3440 type(struct_interpInfo), intent(inout) :: interpInfo
3441 type(struct_gsv) , intent(in) :: stateVector
3442 integer , intent(in) :: headerIndex, kIndex, stepIndex, procIndex
3443 integer , intent(out) :: numGridpt(interpInfo%hco%numSubGrid)
3444
3445 ! Locals:
3446 integer :: depotIndex
3447 integer :: ierr
3448 integer :: latIndex, lonIndex
3449 integer :: subGridIndex
3450 real(4) :: lon_deg_r4, lat_deg_r4
3451 real(4) :: xpos_r4, ypos_r4, xpos2_r4, ypos2_r4
3452
3453 if ( stateVector%hco%grtyp == 'U' ) then
3454 call utl_abort('s2c_setupNearestNeighbor: Yin-Yang grid not supported')
3455 end if
3456
3457 numGridpt(:) = 0
3458
3459 lat_deg_r4 = real(interpInfo%stepProcData(procIndex, stepIndex)%allLat(headerIndex, kIndex) * &
3460 MPC_DEGREES_PER_RADIAN_R8)
3461 lon_deg_r4 = real(interpInfo%stepProcData(procIndex, stepIndex)%allLon(headerIndex, kIndex) * &
3462 MPC_DEGREES_PER_RADIAN_R8)
3463
3464 ierr = gpos_getPositionXY( stateVector%hco%EZscintID, &
3465 xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
3466 lat_deg_r4, lon_deg_r4, subGridIndex )
3467
3468 latIndex = nint(ypos_r4)
3469 lonIndex = nint(xpos_r4)
3470
3471 ! Handle periodicity in longitude
3472 if ( lonIndex == statevector%ni+1 .and. stateVector%hco%grtyp == 'G' ) lonIndex = 1
3473
3474 ! Test bounds
3475 if ( lonIndex < 1 .or. lonIndex > statevector%ni .or. &
3476 latIndex < 1 .or. latIndex > statevector%nj ) then
3477
3478 write(*,*) 's2c_setupNearestNeighbor: observation out of bounds'
3479
3480 else
3481
3482 if ( allocated(interpInfo%interpWeightDepot) ) then
3483
3484 depotIndex = interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex)
3485
3486 interpInfo%interpWeightDepot(depotIndex) = 1.d0
3487 interpInfo%latIndexDepot (depotIndex) = latIndex
3488 interpInfo%lonIndexDepot (depotIndex) = lonIndex
3489
3490 end if
3491
3492 numGridpt(subGridIndex) = 1
3493
3494 end if
3495
3496 end subroutine s2c_setupNearestNeighbor
3497
3498 !--------------------------------------------------------------------------
3499 ! checkColumnStatevectorMatch
3500 !--------------------------------------------------------------------------
3501 subroutine checkColumnStatevectorMatch(column,statevector)
3502 !
3503 !:Purpose: To check column and statevector have identical nk and variables.
3504 !
3505 implicit none
3506
3507 ! Arguments:
3508 type(struct_gsv) , intent(in) :: statevector
3509 type(struct_columnData), intent(in) :: column
3510
3511 ! Locals:
3512 integer :: kIndex
3513
3514 ! check column/statevector have same nk
3515 if ( column%nk /= gsv_getNumK(statevector) ) then
3516 write(*,*) 'checkColumnStatevectorMatch: column%nk, gsv_getNumK(statevector)', column%nk, gsv_getNumK(statevector)
3517 call utl_abort('checkColumnStatevectorMatch: column%nk /= gsv_getNumK(statevector)')
3518 end if
3519
3520 ! loop through k and check varNames are same between column/statevector
3521 do kIndex = 1, column%nk
3522 if (gsv_getVarNameFromK(statevector,kIndex) /= col_getVarNameFromK(column,kIndex)) then
3523 write(*,*) 'checkColumnStatevectorMatch: kIndex, varname in statevector and column: ', kIndex, &
3524 gsv_getVarNameFromK(statevector,kIndex), col_getVarNameFromK(column,kIndex)
3525 call utl_abort('checkColumnStatevectorMatch: varname in column and statevector do not match')
3526 end if
3527 end do
3528
3529 end subroutine checkColumnStatevectorMatch
3530
3531 !--------------------------------------------------------------------------
3532 ! latlonChecks
3533 !--------------------------------------------------------------------------
3534 subroutine latlonChecks( obsSpaceData, hco, headerIndex, rejectOutsideObs, &
3535 latLev_T, lonLev_T, latLev_M, lonLev_M, latLev_S_opt, lonLev_S_opt )
3536 !
3537 !:Purpose: To check if the obs are inside the domain.
3538 !
3539 implicit none
3540
3541 ! Arguments:
3542 type(struct_obs), intent(inout) :: obsSpaceData
3543 type(struct_hco), intent(in) :: hco
3544 integer, intent(in) :: headerIndex
3545 logical, intent(in) :: rejectOutsideObs
3546 real(8), intent(inout) :: latLev_T(:)
3547 real(8), intent(inout) :: lonLev_T(:)
3548 real(8), intent(inout) :: latLev_M(:)
3549 real(8), intent(inout) :: lonLev_M(:)
3550 real(8), optional, intent(inout) :: latLev_S_opt
3551 real(8), optional, intent(inout) :: lonLev_S_opt
3552
3553 ! Locals:
3554 integer :: ierr
3555 integer :: bodyIndex, bodyIndexBeg, bodyIndexEnd, niP1, subGridIndex
3556 integer :: nlev_T, nlev_M
3557 real(4) :: lon_r4, lat_r4, lon_deg_r4, lat_deg_r4
3558 real(4) :: xpos_r4, ypos_r4, xpos2_r4, ypos2_r4
3559 logical :: latlonOutsideGrid, rejectHeader
3560 integer :: gdllfxy
3561
3562 ! Allow for periodicity in Longitude for global Gaussian grid
3563 if ( hco%grtyp == 'G' .or. (hco%grtyp == 'Z' .and. hco%global) ) then
3564 niP1 = hco%ni + 1
3565 else
3566 niP1 = hco%ni
3567 end if
3568
3569 nlev_T = size(latLev_T)
3570 nlev_M = size(latLev_M)
3571
3572 ! check if lat/lon of last thermo level is outside domain.
3573 rejectHeader = .false.
3574 lat_r4 = real(latLev_T(nlev_T),4)
3575 lon_r4 = real(lonLev_T(nlev_T),4)
3576
3577 lat_deg_r4 = lat_r4 * MPC_DEGREES_PER_RADIAN_R8
3578 lon_deg_r4 = lon_r4 * MPC_DEGREES_PER_RADIAN_R8
3579 ierr = gpos_getPositionXY( hco%EZscintID, &
3580 xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
3581 lat_deg_r4, lon_deg_r4, subGridIndex )
3582
3583 latlonOutsideGrid = ( xpos_r4 < 1.0 .or. &
3584 xpos_r4 > real(niP1) .or. &
3585 ypos_r4 < 1.0 .or. &
3586 ypos_r4 > real(hco%nj) )
3587
3588 if ( latlonOutsideGrid .and. rejectOutsideObs ) then
3589 rejectHeader = .true.
3590 end if
3591
3592 ! check if lat/lon of last momentum level is outside domain.
3593 if ( .not. rejectHeader ) then
3594 lat_r4 = real(latLev_M(nlev_M),4)
3595 lon_r4 = real(lonLev_M(nlev_M),4)
3596
3597 lat_deg_r4 = lat_r4 * MPC_DEGREES_PER_RADIAN_R8
3598 lon_deg_r4 = lon_r4 * MPC_DEGREES_PER_RADIAN_R8
3599 ierr = gpos_getPositionXY( hco%EZscintID, &
3600 xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
3601 lat_deg_r4, lon_deg_r4, subGridIndex )
3602
3603 latlonOutsideGrid = ( xpos_r4 < 1.0 .or. &
3604 xpos_r4 > real(niP1) .or. &
3605 ypos_r4 < 1.0 .or. &
3606 ypos_r4 > real(hco%nj) )
3607
3608 if ( latlonOutsideGrid .and. rejectOutsideObs ) then
3609 rejectHeader = .true.
3610 end if
3611 end if
3612
3613 ! check if lat/lon of surface level is outside domain.
3614 if ( present(latLev_S_opt) .and. present(lonLev_S_opt) .and. .not. rejectHeader ) then
3615 lat_r4 = real(latLev_S_opt,4)
3616 lon_r4 = real(lonLev_S_opt,4)
3617
3618 lat_deg_r4 = lat_r4 * MPC_DEGREES_PER_RADIAN_R8
3619 lon_deg_r4 = lon_r4 * MPC_DEGREES_PER_RADIAN_R8
3620 ierr = gpos_getPositionXY( hco%EZscintID, &
3621 xpos_r4, ypos_r4, xpos2_r4, ypos2_r4, &
3622 lat_deg_r4, lon_deg_r4, subGridIndex )
3623
3624 latlonOutsideGrid = ( xpos_r4 < 1.0 .or. &
3625 xpos_r4 > real(niP1) .or. &
3626 ypos_r4 < 1.0 .or. &
3627 ypos_r4 > real(hco%nj) )
3628
3629 if ( latlonOutsideGrid .and. rejectOutsideObs ) then
3630 rejectHeader = .true.
3631 end if
3632 end if
3633
3634 if ( rejectHeader ) then
3635 ! The observation is outside the domain.
3636 ! With a LAM trial field we must discard this observation
3637 write(*,*) 'latlonChecks: Rejecting OBS outside the hco domain, ', headerIndex
3638 write(*,*) ' position : ', lat_deg_r4, lon_deg_r4, ypos_r4, xpos_r4
3639
3640 bodyIndexBeg = obs_headElem_i(obsSpaceData, OBS_RLN, headerIndex)
3641 bodyIndexEnd = obs_headElem_i(obsSpaceData, OBS_NLV, headerIndex) + bodyIndexBeg -1
3642 do bodyIndex = bodyIndexBeg, bodyIndexEnd
3643 call obs_bodySet_i(obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
3644 end do
3645 call obs_headSet_i(obsSpaceData, OBS_ST1, headerIndex, &
3646 ibset( obs_headElem_i(obsSpaceData, OBS_ST1, headerIndex), 05))
3647
3648 ! Assign domain mid-point lat-lon to this header
3649 if ( hco%grtyp == 'Y' ) then
3650 lat_deg_r4 = hco%lat2d_4(hco%ni/2,hco%nj/2)
3651 lon_deg_r4 = hco%lon2d_4(hco%ni/2,hco%nj/2)
3652 else
3653 xpos_r4 = real(hco%ni)/2.0
3654 ypos_r4 = real(hco%nj)/2.0
3655 ierr = gdllfxy(hco%EZscintID, lat_deg_r4, lon_deg_r4, &
3656 xpos_r4, ypos_r4, 1)
3657 end if
3658
3659 lonLev_T(:) = real(lon_deg_r4 * MPC_RADIANS_PER_DEGREE_R4,8)
3660 latLev_T(:) = real(lat_deg_r4 * MPC_RADIANS_PER_DEGREE_R4,8)
3661 lonLev_M(:) = real(lon_deg_r4 * MPC_RADIANS_PER_DEGREE_R4,8)
3662 latLev_M(:) = real(lat_deg_r4 * MPC_RADIANS_PER_DEGREE_R4,8)
3663 if (present(lonLev_S_opt) .and. present(latLev_S_opt)) then
3664 lonLev_S_opt = real(lon_deg_r4 * MPC_RADIANS_PER_DEGREE_R4,8)
3665 latLev_S_opt = real(lat_deg_r4 * MPC_RADIANS_PER_DEGREE_R4,8)
3666 end if
3667
3668 end if
3669
3670 end subroutine latlonChecks
3671
3672 !--------------------------------------------------------------------------
3673 ! getTovsFootprintRadius
3674 !--------------------------------------------------------------------------
3675 function getTovsFootprintRadius(obsSpaceData, headerIndex, beSilent_opt) result(footPrintRadius_r4)
3676 !
3677 !:Purpose: calculate foot-print radius for TOVS observations
3678 !
3679 implicit none
3680
3681 ! Arguments:
3682 type(struct_obs), intent(in) :: obsSpaceData
3683 integer , intent(in) :: headerIndex
3684 logical, optional, intent(in) :: beSilent_opt
3685 ! Result:
3686 real(4) :: footPrintRadius_r4
3687
3688 ! Locals:
3689 integer :: codtyp, sensorIndex
3690 real(8) :: fovAngularDiameter, satHeight, footPrintRadius
3691 character(len=codtyp_name_length) :: instrumName
3692 logical :: beSilent
3693
3694 if ( present(beSilent_opt) ) then
3695 beSilent = beSilent_opt
3696 else
3697 beSilent = .true.
3698 end if
3699
3700 ! get nominal satellite height
3701 sensorIndex = tvs_lsensor(tvs_tovsIndex(headerIndex))
3702 satHeight = tvs_coefs(sensorIndex)%coef%fc_sat_height
3703
3704 ! FOV angular diameter
3705 codtyp = obs_headElem_i( obsSpaceData, OBS_ITY, headerIndex )
3706 instrumName = codtyp_get_name(codtyp)
3707 select case(trim(instrumName))
3708 case('amsua')
3709 fovAngularDiameter = 3.3d0
3710 case('amsub')
3711 fovAngularDiameter = 1.1d0
3712 case('mhs')
3713 fovAngularDiameter = 10.0d0 / 9.0d0
3714 case('airs')
3715 fovAngularDiameter = 1.1d0
3716 case('iasi')
3717 fovAngularDiameter = 14.65d0 / 1000.0d0 * MPC_DEGREES_PER_RADIAN_R8
3718 case('radianceclear')
3719 fovAngularDiameter = 0.125d0
3720 case('ssmis')
3721 fovAngularDiameter = 1.2d0
3722 case('atms')
3723 fovAngularDiameter = 2.2d0
3724 case('cris')
3725 fovAngularDiameter = 14.0d0 / 824.0d0 * MPC_DEGREES_PER_RADIAN_R8
3726 case default
3727 fovAngularDiameter = -1.0d0
3728 end select
3729
3730 if ( fovAngularDiameter < 0.0d0 ) then
3731 footPrintRadius_r4 = bilinearFootprint
3732 else
3733 ! get foot print radius (meter) from angular diameter
3734 footPrintRadius = 0.5d0 * fovAngularDiameter * MPC_RADIANS_PER_DEGREE_R8 * satHeight * 1000
3735 footPrintRadius_r4 = real(footPrintRadius,4)
3736 end if
3737
3738 if ( .not. beSilent ) then
3739 write(*,*) 'getTovsFootprintRadius: sensorIndex=', sensorIndex, &
3740 ',satHeight=', satHeight, ',fovAngularDiameter=', fovAngularDiameter, ',codtyp=', codtyp, &
3741 ',footPrintRadius=', footPrintRadius_r4
3742 end if
3743
3744 end function getTovsFootprintRadius
3745
3746 ! -------------------------------------------------------------
3747 ! s2c_getWeightsAndGridPointIndexes
3748 ! -------------------------------------------------------------
3749 subroutine s2c_getWeightsAndGridPointIndexes(headerIndex, kIndex, stepIndex, procIndex, &
3750 interpWeight, latIndex, lonIndex, gridptCount)
3751 !:Purpose: Returns the weights and grid point indexes for a single observation.
3752 !
3753 !
3754 implicit none
3755
3756 ! Arguments:
3757 integer, intent(in) :: headerIndex
3758 integer, intent(in) :: kIndex
3759 integer, intent(in) :: stepIndex
3760 integer, intent(in) :: procIndex
3761 real(8), intent(out) :: interpWeight(:)
3762 integer, intent(out) :: latIndex(:)
3763 integer, intent(out) :: lonIndex(:)
3764 integer, intent(out) :: gridptCount
3765
3766 ! Locals:
3767 integer :: indexBeg, indexEnd, gridptIndex
3768 integer :: subGridIndex, maxGridpt
3769
3770 call utl_tmg_start(30,'--StateToColumn')
3771
3772 maxGridpt = size( interpWeight )
3773
3774 gridptCount = 0
3775
3776 if ( interpInfo_tlad%stepProcData(procIndex, stepIndex)%allHeaderIndex(headerIndex) /= headerIndex ) then
3777 call utl_abort('s2c_getWeightsAndGridPointIndexes: headerUsedIndex and headerIndex differ.'// &
3778 ' If using multiple time steps in the assimilation window,'// &
3779 ' the code needs to be modified to convert values of headerIndex into headerUsedIndex.')
3780 end if
3781
3782 subGrid_loop: do subGridIndex = 1, interpInfo_tlad%hco%numSubGrid
3783
3784 indexBeg = interpInfo_tlad%stepProcData(procIndex,stepIndex)%depotIndexBeg(subGridIndex, headerIndex, kIndex)
3785 indexEnd = interpInfo_tlad%stepProcData(procIndex,stepIndex)%depotIndexEnd(subGridIndex, headerIndex, kIndex)
3786
3787 if ( indexEnd < IndexBeg ) cycle subGrid_loop
3788
3789 do gridptIndex = indexBeg, indexEnd
3790
3791 gridptCount = gridptCount + 1
3792
3793 if ( gridptCount > maxGridpt ) then
3794 call utl_abort('s2c_getWeightsAndGridPointIndexes: maxGridpt must be increased')
3795 end if
3796
3797 lonIndex(gridptCount) = interpInfo_tlad%lonIndexDepot(gridptIndex)
3798 latIndex(gridptCount) = interpInfo_tlad%latIndexDepot(gridptIndex)
3799 interpWeight(gridptCount) = interpInfo_tlad%interpWeightDepot(gridptIndex)
3800
3801 end do
3802
3803 end do subGrid_loop
3804
3805 call utl_tmg_stop(30)
3806
3807 end subroutine s2c_getWeightsAndGridPointIndexes
3808
3809 ! -------------------------------------------------------------
3810 ! s2c_deallocInterpInfo
3811 ! -------------------------------------------------------------
3812 subroutine s2c_deallocInterpInfo( inputStateVectorType )
3813 !:Purpose: Deallocate interpInfo_nl/tlad object.
3814 !
3815 implicit none
3816
3817 ! Arguments:
3818 character(len=*), intent(in) :: inputStateVectorType
3819
3820 ! Locals:
3821 type(struct_interpInfo), pointer :: interpInfo
3822 integer :: stepIndex, procIndex, numStep
3823
3824 select case( trim(inputStateVectorType) )
3825 case('nl')
3826 interpInfo => interpInfo_nl
3827 case('tlad')
3828 interpInfo => interpInfo_tlad
3829 case default
3830 call utl_abort('s2c_deallocInterpInfo: invalid input argument' // inputStateVectorType)
3831 end select
3832
3833 if ( .not. interpInfo%initialized ) return
3834
3835 write(*,*) 's2c_deallocInterpInfo: deallocating interpInfo for inputStateVectorType=', &
3836 inputStateVectorType
3837
3838 numStep = size(interpInfo%stepProcData,2)
3839
3840 deallocate(interpInfo%interpWeightDepot)
3841 deallocate(interpInfo%latIndexDepot)
3842 deallocate(interpInfo%lonIndexDepot)
3843 do stepIndex = 1, numStep
3844 do procIndex = 1, mmpi_nprocs
3845 deallocate(interpInfo%stepProcData(procIndex,stepIndex)%allLat)
3846 deallocate(interpInfo%stepProcData(procIndex,stepIndex)%allLon)
3847 deallocate(interpInfo%stepProcData(procIndex,stepIndex)%allHeaderIndex)
3848 deallocate(interpInfo%stepProcData(procIndex,stepIndex)%depotIndexBeg)
3849 deallocate(interpInfo%stepProcData(procIndex,stepIndex)%depotIndexEnd)
3850 if ( interpInfo%hco%rotated ) then
3851 deallocate(interpInfo%stepProcData(procIndex,stepIndex)%allLonRot)
3852 deallocate(interpInfo%stepProcData(procIndex,stepIndex)%allLatRot)
3853 end if
3854 end do
3855 end do
3856 deallocate(interpInfo%stepProcData)
3857 deallocate(interpInfo%allNumHeaderUsed)
3858 call oti_deallocate(interpInfo%oti)
3859
3860 interpInfo%initialized = .false.
3861
3862 end subroutine s2c_deallocInterpInfo
3863
3864end module stateToColumn_mod
