ROMS
Loading...
Searching...
No Matches
Functions/Subroutines
rep_matrix.F File Reference
#include "cppdefs.h"
Include dependency graph for rep_matrix.F:

Go to the source code of this file.

Functions/Subroutines

subroutine rep_matrix (ng, model, outloop, ninnloop)
 

Function/Subroutine Documentation

◆ rep_matrix()

subroutine rep_matrix ( integer, intent(in) ng,
integer, intent(in) model,
integer, intent(in) outloop,
integer, intent(in) ninnloop )

Definition at line 3 of file rep_matrix.F.

4!
5!git $Id$
6!================================================== Hernan G. Arango ===
7! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md !
10!=======================================================================
11! !
12! This routine estimates the TRANSPOSE of the inverse of the !
13! stabilized representer matrix using the Lanczos vectors from !
14! the inner-loops. !
15! !
16!=======================================================================
17!
18 USE mod_param
19 USE mod_parallel
20 USE mod_fourdvar
21 USE mod_iounits
22 USE mod_scalars
23
24# ifdef DISTRIBUTE
25!
26 USE distribute_mod, ONLY : mp_bcastf
27# endif
28!
29 implicit none
30!
31! Imported variable declarations
32!
33 integer, intent(in) :: ng, model, outLoop, NinnLoop
34!
35! Local variable declarations.
36!
37 integer :: i, j, iobs, ivec, innLoop
38# ifdef IMPACT_INNER
39 integer :: mloop
40# endif
41# ifdef MINRES
42 integer :: info
43 integer, dimension(NinnLoop) :: ipiv
44# endif
45
46 real(r8) :: zbet
47
48 real(r8), dimension(NinnLoop) :: zu, zgam, zt
49# ifdef MINRES
50 real(r8) :: zsum, zck, zgk
51 real(r8), dimension(NinnLoop,NinnLoop) :: ztriT, zLT, zLTt
52 real(r8), dimension(NinnLoop) :: tau, zwork1, zeref
53# endif
54# ifdef RPCG
55 real(r8), dimension(NinnLoop) :: zfact
56# endif
57# ifdef IMPACT_INNER
58!
59!-----------------------------------------------------------------------
60! Compute observation impact on the weak constraint 4DVAR data
61! assimilation systems for each inner-loop.
62!-----------------------------------------------------------------------
63!
64 master_thread : IF (master) THEN
65
66 DO mloop=1,ninnloop
67
68 DO iobs=1,ndatum(ng)
69 ad_obsval(iobs,mloop)=0.0_r8
70 END DO
71 DO innloop=1,ninnloop
72 zt(innloop)=0.0_r8
73 END DO
74# ifdef RPCG
75 DO innloop=1,ninnloop
76 IF (innloop.eq.1) THEN
77 zfact(innloop)=1.0_r8/cg_qg(1,outloop)
78 ELSE
79 zfact(innloop)=1.0_r8/cg_beta(innloop,outloop)
80 END IF
81 END DO
82# endif
83!
84! Multiply TLmodVal by the tranpose matrix of Lanczos vectors.
85! Note that the factor of 1/SQRT(ObsErr) is added to convert to
86! x-space.
87!
88 DO innloop=1,mloop
89 DO iobs=1,ndatum(ng)
90 IF (obserr(iobs).NE.0.0_r8) THEN
91# ifdef RPCG
92 zt(innloop)=zt(innloop)+ &
93 & obsscale(iobs)* &
94 & zcglwk(iobs,innloop,outloop)* &
95 & tlmodval(iobs)
96# else
97 zt(innloop)=zt(innloop)+ &
98 & obsscale(iobs)* &
99 & zcglwk(iobs,innloop,outloop)* &
100 & tlmodval(iobs)/ &
101 & sqrt(obserr(iobs))
102# endif
103 END IF
104 END DO
105 END DO
106
107# ifdef MINRES
108!
109! Use the minimum residual method as described by Paige and Saunders
110! ("Sparse Indefinite Systems of Linear Equations", 1975, SIAM Journal
111! on Numerical Analysis, 617-619). Specifically we refer to equations
112! 6.10 and 6.11 of this paper.
113!
114! Perform a LQ factorization of the tridiagonal matrix.
115!
116 ztrit=0.0_r8
117 DO i=1,mloop
118 ztrit(i,i)=cg_delta(i,outloop)
119 END DO
120 DO i=1,mloop-1
121 ztrit(i,i+1)=cg_beta(i+1,outloop)
122 END DO
123 DO i=2,mloop
124 ztrit(i,i-1)=cg_beta(i,outloop)
125 END DO
126 CALL sqlq (mloop, ztrit, tau, zwork1)
127!
128! Isolate L=zLT and its tranpose.
129!
130 zlt=0.0_r8
131 zltt=0.0_r8
132 DO i=1,mloop
133 DO j=1,i
134 zlt(i,j)=ztrit(i,j)
135 END DO
136 END DO
137 DO j=1,mloop
138 DO i=1,mloop
139 zltt(i,j)=zlt(j,i)
140 END DO
141 END DO
142!
143! Compute zck.
144!
145 zgk=sqrt(zlt(mloop,mloop)*zlt(mloop,mloop)+ &
146 & cg_beta(mloop+1,outloop)*cg_beta(mloop+1,outloop))
147 zck=zlt(mloop,mloop)/zgk
148!
149! Now form inv(L)*zt - NOTE: we use L not L', so we use the
150! adjoint of the original solver.
151!
152 DO j=1,mloop
153 DO i=j-1,1,-1
154 zt(j)=zt(j)-zt(i)*zltt(i,j)
155 END DO
156 zt(j)=zt(j)/zltt(j,j)
157 END DO
158!
159! Now form zt=D*zt.
160!
161 zt(mloop)=zck*zt(mloop)
162!
163! Now form zt=Q'*zt.
164!
165 DO i=mloop,1,-1
166 DO j=1,mloop
167 zeref(j)=0.0_r8
168 END DO
169 zeref(i)=1.0_r8
170 DO j=i+1,mloop
171 zeref(j)=ztrit(i,j)
172 END DO
173 zsum=0.0_r8
174 DO j=1,mloop
175 zsum=zsum+zt(j)*zeref(j)
176 END DO
177 DO j=1,mloop
178 zt(j)=zt(j)-tau(i)*zsum*zeref(j)
179 END DO
180 END DO
181!
182! Copy the solution zt into zu.
183!
184 DO i=1,mloop
185 zu(i)=zt(i)
186 END DO
187# else
188!
189! Now multiply the result by the inverse tridiagonal matrix.
190!
191 zbet=cg_delta(1,outloop)
192 zu(1)=zt(1)/zbet
193 DO ivec=2,mloop
194 zgam(ivec)=cg_beta(ivec,outloop)/zbet
195 zbet=cg_delta(ivec,outloop)-cg_beta(ivec,outloop)*zgam(ivec)
196 zu(ivec)=(zt(ivec)-cg_beta(ivec,outloop)* &
197 & zu(ivec-1))/zbet
198 END DO
199
200 DO ivec=mloop-1,1,-1
201 zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
202 END DO
203# endif
204!
205! Finally multiply by the matrix of Lanczos vactors.
206# ifndef RPCG
207! Note that the factor of 1/SQRT(ObsErr) is added to covert to
208! x-space.
209# endif
210!
211 DO iobs=1,ndatum(ng)
212 DO innloop=1,mloop
213 IF (obserr(iobs).NE.0.0_r8) THEN
214# ifdef RPCG
215 ad_obsval(iobs,mloop)=ad_obsval(iobs,mloop)+ &
216 & obsscale(iobs)* &
217 & tlmodval_s(iobs,innloop,outloop)* &
218 & zu(innloop)*zfact(innloop)/ &
219 & obserr(iobs)
220# else
221 ad_obsval(iobs,mloop)=ad_obsval(iobs,mloop)+ &
222 & obsscale(iobs)* &
223 & zcglwk(iobs,innloop,outloop)* &
224 & zu(innloop)/ &
225 & sqrt(obserr(iobs))
226# endif
227 END IF
228 END DO
229 END DO
230
231 END DO
232
233 END IF master_thread
234
235# else
236!
237!-----------------------------------------------------------------------
238! Compute observation impact on the weak constraint 4DVAR data
239! assimilation systems.
240!-----------------------------------------------------------------------
241!
242 master_thread : IF (master) THEN
243
244 DO iobs=1,ndatum(ng)
245 ad_obsval(iobs)=0.0_r8
246 END DO
247 DO innloop=1,ninnloop
248 zt(innloop)=0.0_r8
249 END DO
250# ifdef RPCG
251 DO innloop=1,ninnloop
252 IF (innloop.eq.1) THEN
253 zfact(innloop)=1.0_r8/cg_qg(1,outloop)
254 ELSE
255 zfact(innloop)=1.0_r8/cg_beta(innloop,outloop)
256 ENDIF
257 END DO
258# endif
259!
260! Multiply TLmodVal by the tranpose matrix of Lanczos vectors.
261! Note that the factor of 1/SQRT(ObsErr) is added to convert to
262! x-space.
263!
264 DO innloop=1,ninnloop
265 DO iobs=1,ndatum(ng)
266 IF (obserr(iobs).NE.0.0_r8) THEN
267# ifdef RPCG
268 zt(innloop)=zt(innloop)+obsscale(iobs)* &
269 & zcglwk(iobs,innloop,outloop)*tlmodval(iobs)
270# else
271 zt(innloop)=zt(innloop)+obsscale(iobs)* &
272 & zcglwk(iobs,innloop,outloop)*tlmodval(iobs)/ &
273 & sqrt(obserr(iobs))
274# endif
275 END IF
276 END DO
277 END DO
278
279# ifdef MINRES
280!
281! Use the minimum residual method as described by Paige and Saunders
282! ("Sparse Indefinite Systems of Linear Equations", 1975, SIAM Journal
283! on Numerical Analysis, 617-619). Specifically we refer to equations
284! 6.10 and 6.11 of this paper.
285!
286! Perform a LQ factorization of the tridiagonal matrix.
287!
288 ztrit=0.0_r8
289 DO i=1,ninnloop
290 ztrit(i,i)=cg_delta(i,outloop)
291 END DO
292 DO i=1,ninnloop-1
293 ztrit(i,i+1)=cg_beta(i+1,outloop)
294 END DO
295 DO i=2,ninnloop
296 ztrit(i,i-1)=cg_beta(i,outloop)
297 END DO
298 CALL sqlq (ninnloop, ztrit, tau, zwork1)
299!
300! Isolate L=zLT and its tranpose.
301!
302 zlt=0.0_r8
303 zltt=0.0_r8
304 DO i=1,ninnloop
305 DO j=1,i
306 zlt(i,j)=ztrit(i,j)
307 END DO
308 END DO
309 DO j=1,ninnloop
310 DO i=1,ninnloop
311 zltt(i,j)=zlt(j,i)
312 END DO
313 END DO
314!
315! Compute zck.
316!
317 zgk=sqrt(zlt(ninnloop,ninnloop)*zlt(ninnloop,ninnloop)+ &
318 & cg_beta(ninnloop+1,outloop)*cg_beta(ninnloop+1,outloop))
319 zck=zlt(ninnloop,ninnloop)/zgk
320!
321! Now form inv(L)*zt - NOTE: we use L not L', so we use the
322! adjoint of the original solver.
323!
324 DO j=1,ninnloop
325 DO i=j-1,1,-1
326 zt(j)=zt(j)-zt(i)*zltt(i,j)
327 END DO
328 zt(j)=zt(j)/zltt(j,j)
329 END DO
330!
331! Now form zt=D*zt.
332!
333 zt(ninnloop)=zck*zt(ninnloop)
334!
335! Now form zt=Q'*zt.
336!
337 DO i=ninnloop,1,-1
338 DO j=1,ninnloop
339 zeref(j)=0.0_r8
340 END DO
341 zeref(i)=1.0_r8
342 DO j=i+1,ninnloop
343 zeref(j)=ztrit(i,j)
344 END DO
345 zsum=0.0_r8
346 DO j=1,ninnloop
347 zsum=zsum+zt(j)*zeref(j)
348 END DO
349 DO j=1,ninnloop
350 zt(j)=zt(j)-tau(i)*zsum*zeref(j)
351 END DO
352 END DO
353!
354! Copy the solution zt into zu.
355!
356 DO i=1,ninnloop
357 zu(i)=zt(i)
358 END DO
359# else
360!
361! Now multiply the result by the inverse tridiagonal matrix.
362!
363 zbet=cg_delta(1,outloop)
364 zu(1)=zt(1)/zbet
365 DO ivec=2,ninnloop
366 zgam(ivec)=cg_beta(ivec,outloop)/zbet
367 zbet=cg_delta(ivec,outloop)-cg_beta(ivec,outloop)*zgam(ivec)
368 zu(ivec)=(zt(ivec)-cg_beta(ivec,outloop)* &
369 & zu(ivec-1))/zbet
370 END DO
371
372 DO ivec=ninnloop-1,1,-1
373 zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
374 END DO
375# endif
376!
377! Finally multiply by the matrix of Lanczos vactors.
378# ifndef RPCG
379! Note that the factor of 1/SQRT(ObsErr) is added to covert to
380! x-space.
381# endif
382!
383 DO iobs=1,ndatum(ng)
384 DO innloop=1,ninnloop
385 IF (obserr(iobs).NE.0.0_r8) THEN
386# ifdef RPCG
387 ad_obsval(iobs)=ad_obsval(iobs)+ &
388 & obsscale(iobs)* &
389 & tlmodval_s(iobs,innloop,outloop)* &
390 & zu(innloop)*zfact(innloop)/ &
391 & obserr(iobs)
392# else
393 ad_obsval(iobs)=ad_obsval(iobs)+ &
394 & obsscale(iobs)* &
395 & zcglwk(iobs,innloop,outloop)* &
396 & zu(innloop)/ &
397 & sqrt(obserr(iobs))
398# endif
399 END IF
400 END DO
401 END DO
402
403 END IF master_thread
404
405# endif
406# ifdef DISTRIBUTE
407!
408!-----------------------------------------------------------------------
409! Broadcast new solution to other nodes.
410!-----------------------------------------------------------------------
411!
412 CALL mp_bcastf (ng, model, ad_obsval)
413# endif
414!
415 RETURN
distribute_mod::mp_bcastf
Definition distribute.F:75
distribute_mod
Definition distribute.F:3
mod_fourdvar
Definition mod_fourdvar.F:2
mod_fourdvar::cg_beta
real(dp), dimension(:,:), allocatable cg_beta
Definition mod_fourdvar.F:654
mod_fourdvar::ndatum
integer, dimension(:), allocatable ndatum
Definition mod_fourdvar.F:409
mod_fourdvar::cg_qg
real(dp), dimension(:,:), allocatable cg_qg
Definition mod_fourdvar.F:706
mod_fourdvar::ad_obsval
real(r8), dimension(:,:), allocatable ad_obsval
Definition mod_fourdvar.F:322
mod_fourdvar::obsscale
real(r8), dimension(:), allocatable obsscale
Definition mod_fourdvar.F:227
mod_fourdvar::obserr
real(r8), dimension(:), allocatable obserr
Definition mod_fourdvar.F:225
mod_fourdvar::zcglwk
real(r8), dimension(:,:,:), allocatable zcglwk
Definition mod_fourdvar.F:288
mod_fourdvar::cg_delta
real(dp), dimension(:,:), allocatable cg_delta
Definition mod_fourdvar.F:691
mod_iounits
Definition mod_iounits.F:2
mod_parallel
Definition mod_parallel.F:2
mod_parallel::master
logical master
Definition mod_parallel.F:40
mod_param
Definition mod_param.F:2
mod_scalars
Definition mod_scalars.F:2

References mod_fourdvar::ad_obsval, mod_fourdvar::cg_beta, mod_fourdvar::cg_delta, mod_fourdvar::cg_qg, mod_parallel::master, mod_fourdvar::ndatum, mod_fourdvar::obserr, mod_fourdvar::obsscale, and mod_fourdvar::zcglwk.

Referenced by roms_kernel_mod::roms_run().

Here is the caller graph for this function: