rMatrix_mod sourceΒΆ

  1module rMatrix_mod
  2  ! MODULE rMatrix_mod (prefix='rmat' category='2. B and R matrices')
  3  !
  4  !:Purpose:  Module to handle non-diagonal observation-error covariance
  5  !           matrices for assimilation of radiances.
  6  !
  7  use rttovInterfaces_mod
  8  use midasMpi_mod
  9  use rttov_const, only  : errorstatus_success
 10  use utilities_mod
 11  use obsSpaceData_mod
 12  use tovsNL_mod
 13  implicit none
 14  private
 15  save
 16
 17  ! public variables
 18  public :: rmat_lnondiagr
 19  ! public subroutines
 20  public :: rmat_init,rmat_cleanup,rmat_readCMatrix,rmat_RsqrtInverseOneObs, rmat_RsqrtInverseAllObs
 21
 22 type rmat_matrix
 23    real(8), pointer     :: Rmat(:,:)=>null()
 24    integer, pointer     :: listChans(:)=>null()
 25    integer              :: nchans=0
 26 END type rmat_matrix
 27
 28  type(rmat_matrix),target,allocatable  :: Rcorr_inst(:) ! non diagonal Correlation matrices for each instrument
 29  type(rmat_matrix),target,allocatable  :: R_tovs(:) ! non diagonal R matrices used for the assimilation of all radiances
 30
 31  ! namelist variable
 32  logical :: rmat_lnondiagr ! choose to use non-diagonal R matrix (i.e. non-zero correlations)
 33
 34  contains
 35
 36  subroutine rmat_init(nsensors,nobtovs)
 37   
 38    implicit none
 39
 40    ! Arguments:
 41    integer, intent(in) :: nsensors
 42    integer, intent(in) :: nobtovs
 43
 44    ! Locals:
 45    integer :: nulnam,ierr
 46    integer, external:: fnom,fclos
 47    namelist /NAMRMAT/rmat_lnonDiagR
 48
 49    ! Default value for parameter rmat_lnondiagr, don't use interchannel correlation by default
 50    rmat_lnonDiagR = .false.
 51
 52    ! Read the parameters from NAMRMAT
 53    nulnam = 0
 54    ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
 55    read(nulnam,nml=namrmat, iostat=ierr)
 56    if (ierr /= 0) call utl_abort('rmat_init: Error reading namelist')
 57    if (mmpi_myid == 0) write(*,nml=namrmat)
 58    ierr = fclos(nulnam)
 59    if (rmat_lnonDiagR) then
 60      allocate(Rcorr_inst(nsensors))
 61      allocate(R_tovs(nobtovs))
 62    end if
 63
 64  end subroutine rmat_init
 65
 66
 67  subroutine rmat_cleanup()
 68    implicit none
 69
 70    if (rmat_lnondiagr) then
 71      deallocate(Rcorr_inst)
 72      deallocate(R_tovs)
 73    end if
 74
 75  end subroutine rmat_cleanup
 76
 77
 78  subroutine rmat_readCMatrix(instrument, sensor_id, ichan )
 79    implicit none
 80    
 81    ! Arguments:
 82    integer, intent(in) :: instrument(3)
 83    integer, intent(in) :: sensor_id
 84    integer, intent(in) :: ichan(:)
 85
 86    ! Locals:
 87    character (len=64) :: filename
 88    integer :: err
 89
 90    call rttov_coeffname (err, instrument, coeffname=filename, filetype="Cmat")
 91    
 92    if (err == errorstatus_success) then
 93       filename = trim(filename) // ".dat"
 94       call rmat_readCMatrixByFileName(filename,Rcorr_inst(sensor_id), ichan )
 95    else
 96      write(*,*) "Unknown instrument ",instrument(:)
 97      call utl_abort("rmat_read_C_matrix")
 98    end if
 99    
100  end subroutine rmat_readCMatrix
101
102
103  subroutine rmat_readCMatrixByFileName(infile,C,chanList_opt)
104    implicit none
105
106    ! Arguments:
107    character (len=*), intent(in)    :: infile          ! name of input file
108    type(rmat_matrix), intent(inout) :: C               ! correlation matrix structure
109    integer, optional, intent(in)    :: chanList_opt(:) ! list of requested channels (if missing will read all file content)
110
111    ! Locals:
112    integer :: i,j,iu,ierr,count,ich,nchn,nch
113    integer, external :: fnom,fclos
114    real(8) :: x
115    integer, allocatable :: index(:)
116
117    nchn = -1
118    if (present(chanList_opt)) then
119      nchn=size(chanList_opt)
120    end if
121
122    iu = 0
123    ierr = fnom(iu,trim(infile),'FTN+SEQ+R/O',0)
124    if (ierr /= 0) then
125      write(*,*) "Cannot open "//trim(infile)
126      call utl_abort("rmat_readCMatrixByFileName")
127    end if
128
129    write(*,*) "rmat_readCMatrixByFileName: Reading "//trim(infile)
130    
131    read(iu,*) nch
132    if (nchn == -1) then
133      nchn = nch
134    else
135      if(nchn > nch) then
136        write(*,*) "Not enough channels in "//trim(infile)
137        write(*,*) nchn,nch
138        call utl_abort("rmat_readCMatrixByFileName")
139      end if
140    end if
141    allocate(index(nch))
142    
143    C%nchans = nchn
144    allocate(C%Rmat(nchn,nchn))
145    allocate(C%listChans(nchn))
146    C%Rmat = 0.d0
147    do i = 1,nchn
148      C%Rmat(i,i) = 1.d0
149    end do
150    count = 0
151    index = -1
152    do i=1,nch
153      read(iu,*) ich
154      if (present(chanList_opt)) then
155        bj:do j=1,nchn
156          if (ich == chanList_opt(j)) then
157            count = count + 1
158            index(i) = j
159            C%listChans(count) = ich
160            exit bj
161          end if
162        end do bj
163      else
164        index(i) = i
165        C%listChans(i) = ich
166        count = count + 1
167      end if
168    end do
169    if (count /= nchn) then
170      write(*,*) "Warning: Missing information in "//trim(infile)
171      do j=1,nchn
172        write(*,*) j, chanList_opt(j) 
173      end do
174      write(*,*) "Not important if there is no observation of this family"
175    end if
176
177    do
178      read(iu,*,iostat=ierr) i,j,x
179      if (ierr /= 0) exit
180      if (index(i) /= -1 .and. index(j) /= -1) then
181        C%Rmat(index(i),index(j)) = x
182        C%Rmat(index(j),index(i)) = x
183      end if
184    end do
185
186    ierr= fclos(iu)
187    deallocate(index)
188
189  end subroutine rmat_readCMatrixByFileName
190
191
192  subroutine rmat_RsqrtInverseOneObs(sensor_id,nsubset,obsIn,obsOut,list_sub,list_oer,indexTovs)
193    !
194    !:Purpose: Apply the operator R**-1/2 to obsIn
195    !           result in obsOut for the subset of channels specified
196    !           in list_sub
197    !
198    implicit none
199
200    ! Arguments:
201    integer, intent(in)  :: sensor_id
202    integer, intent(in)  :: nsubset
203    integer, intent(in)  :: list_sub(nsubset)
204    real(8), intent(in)  :: list_oer(nsubset)
205    real(8), intent(in)  :: obsIn(nsubset)
206    real(8), intent(out) :: obsOut(nsubset)
207    integer, intent(in)  :: indexTovs
208
209    ! Locals:
210    real (8) :: Rsub(nsubset,nsubset), alpha, beta, product 
211    integer :: index(nsubset)
212    integer :: i,j
213
214    if (R_tovs(indexTovs)%nchans == 0) then
215
216      if (sensor_id <= 0 .or. sensor_id > size(Rcorr_inst)) then
217        write(*,*) "invalid sensor_id",sensor_id,size(Rcorr_inst)
218        call utl_abort('rmat_RsqrtInverseOneObs')
219      end if
220
221      index = -1
222      do i=1,nsubset
223        bj: do j = 1, Rcorr_inst(sensor_id)%nchans
224          if (list_sub(i) == Rcorr_inst(sensor_id)%listChans(j)) then
225            index(i) = j
226            exit bj 
227          end if
228        end do bj
229      end do
230      if (any(index == -1)) then
231        write(*,*) "Missing information for some channel !"
232        write(*,*) list_sub(:)
233        write(*,*) index(:)
234        call utl_abort('rmat_RsqrtInverseOneObs')
235      end if
236      R_tovs(indexTovs)%nchans = nsubset
237      allocate(R_tovs(indexTovs)%listChans(nsubset))
238      R_tovs(indexTovs)%listChans(1:nsubset) = list_sub(1:nsubset)
239      do j=1,nsubset
240        do i=1,nsubset
241          product = list_oer(i) * list_oer(j)
242          Rsub(i,j) = product * Rcorr_inst(sensor_id)%Rmat(index(i),index(j))
243        end do
244      end do
245      ! Calculation of R**-1/2
246      call utl_tmg_start(20,'----RmatMatSqrt')
247      call utl_matSqrt(Rsub,nsubset,-1.d0,.false.)
248      call utl_tmg_stop(20)
249      allocate(R_tovs(indexTovs)%Rmat(nsubset,nsubset))
250      do j=1,nsubset
251        do i=1,nsubset
252          R_tovs(indexTovs)%Rmat(i,j) = Rsub(i,j)
253        end do
254      end do
255    end if
256
257    call utl_tmg_start(21,'----RmatMatMult')
258    alpha = 1.d0
259    beta = 0.d0
260    obsOut = 0.d0
261    ! Optimized symetric matrix vector product from Lapack
262    call dsymv("L", nsubset, alpha, R_tovs(indexTovs)%Rmat, nsubset,obsIn, 1, beta, obsOut, 1)
263    call utl_tmg_stop(21)
264
265  end subroutine rmat_RsqrtInverseOneObs
266
267
268  !--------------------------------------------------------------------------
269  ! rmat_RsqrtInverseAllObs
270  !--------------------------------------------------------------------------
271  subroutine rmat_RsqrtInverseAllObs( obsSpaceData, elem_dest_i, elem_src_i )
272    !
273    !:Purpose: To apply observation-error variances to ROBDATA8(k_src,*) and to
274    !          store it in the elem_src_s of obsspacedata
275    !
276    implicit none
277
278    ! Arguments:
279    type(struct_obs), intent(inout) :: obsspacedata
280    integer,          intent(in)    :: elem_dest_i ! destination index
281    integer,          intent(in)    :: elem_src_i  ! source index
282
283    ! Locals:
284    integer :: bodyIndex, headerIndex
285    integer :: idata, idatend, idatyp, count, channelNumber, channelIndex
286    real(8) :: obsIn( tvs_maxChannelNumber ), obsOut( tvs_maxChannelNumber )
287    integer :: list_chan( tvs_maxChannelNumber )
288    real(8) :: list_OER( tvs_maxChannelNumber )
289
290    ! NOTE I tried using openMP on this loop, but it increased the cost from 4sec to 80sec!!!
291    do headerIndex =1, obs_numHeader(obsspacedata)
292
293      idata   = obs_headElem_i( obsspacedata, OBS_RLN, headerIndex )
294      idatend = obs_headElem_i( obsspacedata, OBS_NLV, headerIndex ) + idata - 1
295      idatyp  = obs_headElem_i( obsspacedata, OBS_ITY, headerIndex )
296
297      if ( tvs_isIdBurpTovs(idatyp) .and. rmat_lnondiagr ) then
298
299        count = 0
300
301        do bodyIndex = idata, idatend
302
303          if (obs_bodyElem_i( obsspacedata, OBS_ASS, bodyIndex ) == obs_assimilated ) then
304            call tvs_getChannelNumIndexFromPPP( obsSpaceData, headerIndex, bodyIndex, &
305                                                channelNumber, channelIndex )
306            count = count + 1
307            list_chan( count ) = channelNumber
308            list_OER( count ) = obs_bodyElem_r( obsspacedata, OBS_OER, bodyIndex )
309            obsIn( count ) = obs_bodyElem_r( obsspacedata, elem_src_i, bodyIndex )
310          end if
311        
312        end do
313
314        if ( count > 0 .and. tvs_tovsIndex( headerIndex ) > 0 ) then
315
316          call rmat_RsqrtInverseOneObs( tvs_lsensor( tvs_tovsIndex( headerIndex )), count, obsIn(1:count), obsOut(1:count), list_chan(1:count), list_OER(1:count), tvs_tovsIndex(headerIndex) )
317
318          count = 0
319          do bodyIndex = idata, idatend
320            if ( obs_bodyElem_i( obsspacedata, OBS_ASS, bodyIndex ) == obs_assimilated) then
321              count = count + 1
322              call obs_bodySet_r(obsspacedata, elem_dest_i, bodyIndex,obsOut(count))
323            end if
324          end do
325
326        else
327
328          do bodyIndex = idata, idatend
329            call obs_bodySet_r(obsspacedata, elem_dest_i, bodyIndex, 0.d0)
330          end do
331
332        end if 
333
334      else
335       
336        do bodyIndex = idata, idatend
337          if (obs_bodyElem_i( obsspacedata, OBS_ASS, bodyIndex ) == obs_assimilated) then
338            call obs_bodySet_r( obsspacedata, elem_dest_i, bodyIndex, &
339                 obs_bodyElem_r( obsspacedata, elem_src_i, bodyIndex) / obs_bodyElem_r( obsspacedata, OBS_OER, bodyIndex ))
340          end if
341        end do
342
343      end if ! is it a radiance in non diagonal R mode ?
344
345    end do !loop on header
346 
347  end subroutine rmat_RsqrtInverseAllObs
348
349end module rMatrix_mod