1module var1D_mod
2 ! MODULE var1D_mod (prefix='var1D' category='4. Data Object transformations')
3 !
4 !:Purpose: contains all 1Dvar-related methods.
5 !
6 use columnData_mod
7 use gridStatevector_mod
8 use horizontalCoord_mod
9 use midasMpi_mod
10 use obsSpaceData_mod
11 use timeCoord_mod
12 use verticalCoord_mod
13 use codeprecision_mod
14 use mathphysconstants_mod
15
16 implicit none
17 save
18 private
19
20 ! public procedures
21 public :: var1D_Setup, var1D_Finalize
22 public :: var1D_transferColumnToYGrid
23 ! public variables
24 public :: var1D_validHeaderIndex, var1D_validHeaderCount
25
26 logical :: initialized = .false.
27 integer, external :: get_max_rss
28 integer, allocatable :: var1D_validHeaderIndex(:) ! pointeur vers les colonnes assimilables pour minimiser la taille du vecteur de controle
29 integer :: var1D_validHeaderCount ! taille effective de var1D_validHeaderIndex
30
31contains
32
33 !--------------------------------------------------------------------------
34 ! var1D_setup
35 !--------------------------------------------------------------------------
36 subroutine var1D_setup(obsSpaceData)
37 !
38 !:Purpose: to setup var1D module
39 !
40 implicit none
41
42 ! Arguments:
43 type (struct_obs), intent(in) :: obsSpaceData
44
45 ! Locals:
46 integer :: countGood, headerIndex
47 integer :: bodyStart, bodyEnd, bodyIndex
48
49 if(mmpi_myid == 0) write(*,*) 'var1D_setup: Starting'
50 if(mmpi_myid == 0) write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
51
52 !we want to count how many obs are really assimilable to minimize controlvector size
53 var1D_validHeaderCount = 0
54 allocate( var1D_validHeaderIndex(obs_numHeader(obsSpaceData)) )
55 do headerIndex = 1, obs_numHeader(obsSpaceData)
56 bodyStart = obs_headElem_i(obsSpaceData, OBS_RLN, headerIndex)
57 bodyEnd = obs_headElem_i(obsSpaceData, OBS_NLV, headerIndex) + bodyStart - 1
58 countGood = 0
59 do bodyIndex = bodyStart, bodyEnd
60 if (obs_bodyElem_i(obsSpaceData, OBS_ASS, bodyIndex) == obs_assimilated) countGood = countGood + 1
61 end do
62 if (countGood > 0) then
63 var1D_validHeaderCount = var1D_validHeaderCount + 1
64 var1D_validHeaderIndex(var1D_validHeaderCount) = headerIndex
65 if (var1D_validHeaderCount == 1) write(*,*) 'first OBS', headerIndex
66 end if
67 end do
68 initialized = .true.
69 write(*,*) 'var1D_setup: var1D_validHeaderCount, obs_numHeader(obsSpaceData)', var1D_validHeaderCount, obs_numHeader(obsSpaceData)
70
71 end subroutine var1D_setup
72
73 !--------------------------------------------------------------------------
74 ! var1D_Finalize
75 !--------------------------------------------------------------------------
76 subroutine var1D_Finalize()
77 !
78 !:Purpose: to deallocate memory used by internal module structures
79 !
80 implicit none
81
82 if (initialized) then
83 deallocate( var1D_validHeaderIndex )
84 end if
85
86 end subroutine var1D_Finalize
87
88 !--------------------------------------------------------------------------
89 ! var1D_transferColumnToYGrid
90 !--------------------------------------------------------------------------
91 subroutine var1D_transferColumnToYGrid( stateVector, obsSpaceData, column, varList)
92 !
93 !:Purpose: to transfer content of a columndata object to a statevector object
94 ! without interpolation (to be used in 1DVar mode to write increments on Y grid).
95 !
96 implicit none
97
98 ! Arguments:
99 type(struct_gsv), intent(inout) :: stateVector
100 type(struct_obs), intent(in) :: obsSpaceData
101 type(struct_columnData), intent(inout) :: column
102 character(len=4), intent(in) :: varList(:)
103
104 ! Locals:
105 type(struct_hco), pointer :: hco_yGrid
106 integer :: varIndex, globalObsIndex, obsIndex, taskIndex, headerIndex
107 integer, allocatable :: var1D_validHeaderCountAllTasks(:), obsOffset(:)
108 real(8), pointer :: myColumn(:), myField(:,:,:)
109 real(8), allocatable, target :: dummy(:)
110 integer :: var1D_validHeaderCountMpiGlobal, var1D_validHeaderCountMax, ierr, status
111 real(8) :: lat, lon
112 integer :: varDim, tag
113
114 call rpn_comm_barrier("GRID",ierr)
115 allocate( obsOffset(0:mmpi_nprocs-1) )
116 if (mmpi_myid ==0) then
117 allocate( var1D_validHeaderCountAllTasks(mmpi_nprocs) )
118 else
119 allocate(var1D_validHeaderCountAllTasks(1))
120 end if
121
122 call rpn_comm_gather(var1D_validHeaderCount , 1, 'MPI_INTEGER', var1D_validHeaderCountAllTasks, 1,'MPI_INTEGER', 0, "GRID", ierr )
123 if (mmpi_myId ==0) then
124 var1D_validHeaderCountMpiGlobal = sum( var1D_validHeaderCountAllTasks(:) )
125 var1D_validHeaderCountMax = maxval( var1D_validHeaderCountAllTasks(:) )
126 obsOffset(0) = 0
127 do taskIndex = 1, mmpi_nprocs - 1
128 obsOffset(taskIndex) = obsOffset(taskIndex - 1) + var1D_validHeaderCountAllTasks(taskIndex)
129 end do
130 write(*,*) 'var1D_transferColumnToYGrid: obsOffset: ', obsOffset(:)
131 end if
132 call rpn_comm_bcast( obsOffset, mmpi_nprocs, 'MPI_INTEGER', 0, "GRID",ierr )
133 call rpn_comm_bcast( var1D_validHeaderCountMax, 1, 'MPI_INTEGER', 0, "GRID",ierr )
134
135 call hco_setupYgrid(hco_Ygrid, 1, var1D_validHeaderCountMpiGlobal)
136 if (mmpi_myId ==0) then
137 call gsv_allocate(stateVector, numstep=tim_nstepobsinc, hco_ptr=hco_Ygrid, vco_ptr=column%vco, &
138 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.false., &
139 dataKind_opt=pre_incrReal, allocHeight_opt=.false., allocPressure_opt=.false., &
140 besilent_opt=.false.)
141 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
142 end if
143
144 write(*,*) 'var1D_transferColumnToYGrid: start of lat-lon dissemination'
145 do obsIndex = 1, var1D_validHeaderCountMax
146 if (obsIndex <= var1D_validHeaderCount ) then
147 headerIndex = var1D_validHeaderIndex(obsIndex)
148 lat = obs_headElem_r(obsSpaceData, OBS_LAT, headerIndex)
149 lon = obs_headElem_r(obsSpaceData, OBS_LON, headerIndex)
150 else
151 lat = MPC_missingValue_R8
152 lon = MPC_missingValue_R8
153 end if
154 if (mmpi_myId == 0) then
155 if ( obsIndex <= var1D_validHeaderCount ) then
156 hco_yGrid%lat2d_4(1, obsIndex) = lat
157 hco_yGrid%lon2d_4(1, obsIndex) = lon
158 end if
159 else
160 tag = 2 * mmpi_myID
161 call rpn_comm_send( lat, 1, 'mpi_real8', 0, tag, 'GRID', ierr )
162 call rpn_comm_send( lon, 1, 'mpi_real8', 0, tag + 1, 'GRID', ierr )
163 end if
164
165 if (mmpi_myId == 0) then
166 do taskIndex = 1, mmpi_nprocs - 1
167 tag = 2 * taskIndex
168 call rpn_comm_recv( lat, 1, 'mpi_real8', taskIndex, tag, 'GRID', status, ierr )
169 call rpn_comm_recv( lon, 1, 'mpi_real8', taskIndex, tag+1, 'GRID', status, ierr )
170 if (lat /= MPC_missingValue_R8 .and. lon /= MPC_missingValue_R8) then
171 globalObsIndex = obsIndex + obsOffset(taskIndex)
172 hco_yGrid%lat2d_4(1, globalObsIndex) = lat
173 hco_yGrid%lon2d_4(1, globalObsIndex) = lon
174 end if
175 end do
176 end if
177 call rpn_comm_barrier("GRID",ierr)
178 end do
179
180 call rpn_comm_barrier("GRID",ierr)
181 write(*,*) 'var1D_transferColumnToYGrid: end of lat-lon dissemination'
182
183 do varIndex = 1, size(varList)
184 write(*,*) 'var1D_transferColumnToYGrid: start of dissemination for ', varList(varIndex)
185 if (mmpi_myId == 0 ) then
186 call gsv_getField(stateVector, myField, varName_opt=varList(varIndex), stepIndex_opt=1)
187 varDim = gsv_getNumLevFromVarName(stateVector, varList(varIndex))
188 end if
189 call rpn_comm_bcast(varDim, 1, 'MPI_INTEGER', 0, 'GRID', ierr)
190 allocate(dummy(varDim))
191 dummy(:) = MPC_missingValue_R8
192 do obsIndex = 1, var1D_validHeaderCountMax
193 if (obsIndex <= var1D_validHeaderCount ) then
194 headerIndex = var1D_validHeaderIndex(obsIndex)
195 myColumn => col_getColumn(column, headerIndex, varName_opt=varList(varIndex))
196 else
197 myColumn => dummy
198 end if
199 if (mmpi_myId == 0) then
200 if ( obsIndex <= var1D_validHeaderCount ) then
201 myField(1, obsIndex, :) = myColumn(:)
202 end if
203 else
204 tag = mmpi_myId
205 call rpn_comm_send(myColumn , varDim, 'mpi_real8', 0, tag, 'GRID', ierr )
206 end if
207
208 if (mmpi_myId == 0) then
209 do taskIndex = 1, mmpi_nprocs - 1
210 tag = taskIndex
211 call rpn_comm_recv(myColumn, varDim, 'mpi_real8', taskIndex, tag, 'GRID', status, ierr )
212 if (all( myColumn /= MPC_missingValue_R8)) then
213 globalObsIndex = obsIndex + obsOffset(taskIndex)
214 myField(1, globalObsIndex, :) = myColumn(:)
215 end if
216 end do
217 end if
218 end do
219
220 write(*,*) 'var1D_transferColumnToYGrid: end of dissemination for ', varList(varIndex)
221 deallocate(dummy)
222
223 end do
224
225 call rpn_comm_barrier("GRID", ierr)
226 deallocate( obsOffset )
227 deallocate( var1D_validHeaderCountAllTasks )
228
229 end subroutine var1D_transferColumnToYGrid
230
231end module var1D_mod
