1program midas_adjointTest
2 !
3 !:Purpose: Main program for adjoint test applications, i.e., testing if the identity
4 ! :math:`<x,L(y)> = <L^T(x),y>` is respected.
5 !
6 !:Algorithm: Initialize a random grid-point space statevector :math:`x` and a random
7 ! control vector :math:`y`. Compute :math:`L(y)` by the action of a linear
8 ! operator and compute the inner product :math:`<x,L(y)>`. Then compute
9 ! :math:`L^T(x)` through the action of the corresponding adjoint operator
10 ! and the inner product :math:`<L^T(x),y>`. Verify if both inner products
11 ! are equal.
12 !
13 ! --
14 !
15 !:File I/O: The required input files vary according to the application (see options below).
16 !
17 !================================================= ==============================================================
18 ! Input and Output Files (square root covariances) Description of file
19 !================================================= ==============================================================
20 ! ``flnml`` In - Main namelist file with parameters user may modify
21 ! ``trlm_01`` In - Background state (a.k.a. trial) file
22 ! ``analysisgrid`` In - File defining grid for computing the analysis increment
23 ! ``bgcov`` In - Static (i.e. NMC) B matrix file for NWP fields
24 ! ``ensemble/$YYYYMMDDHH_006_$NNNN`` In - Ensemble member files defining ensemble B matrix
25 ! ``innerProd.txt`` Out - Results of inner products and their difference
26 !================================================= ==============================================================
27 !
28 !:Synopsis: There are several flavors of tests that can be run. The choice is configured
29 ! in the namelist block ``NAMADT`` with the variable ``test`` (see below).
30 ! All tests follow this general synopsis:
31 !
32 ! - initialize with gaussian noise (``rng_gaussian()``) a statevector
33 ! :math:`x` and a control vector :math:`y`.
34 !
35 ! - apply the corresponding square root operator on :math:`y` to obtain :math:`Ly`
36 !
37 ! - compute the inner product :math:`<x ,Ly>`
38 !
39 ! - apply the corresponding square root adjoint operator on :math:`x` to
40 ! obtain :math:`L^T(x)`
41 !
42 ! - compute the inner product :math:`<L^Tx,y>`
43 !
44 ! - print :math:`<x,L(y)>`, :math:`<L^T(x),y>` and their relative difference
45 !
46 !:Options: `List of namelist blocks <../namelists_in_each_program.html#adjointtest>`_
47 ! that can affect the ``adjointTest`` program.
48 !
49 ! * As described in the algorithm section, four different tests are implemented.
50 ! The test that is conducted is chosen through the namelist variable
51 ! `&NAMADT Test` that can be either
52 !
53 ! * **'Bhi'** : test the adjoint of the square root of the homogeneous and
54 ! isotropic covariance matrix. The namelist block `&NAMBHI` will
55 ! configure the covariance matrix properties.
56 !
57 ! * **'Bens'** : test the adjoint of the square root of the ensemble covariance
58 ! matrix. The namelist block `&NAMBEN` will configure the covariance
59 ! matrix properties.
60 !
61 ! * **'loc'** : test the adjoint of the square root of localized covariance
62 ! matrix.
63 !
64 ! * **'advEns'** : test the adjoint of the advection operator on ensemble
65 !
66 ! * **'advGSV'** : test the adjoint of the advection operator on a single
67 ! statevector
68 !
69 use version_mod
70 use codePrecision_mod
71 use ramDisk_mod
72 use utilities_mod
73 use midasMpi_mod
74 use mathPhysConstants_mod
75 use horizontalCoord_mod
76 use verticalCoord_mod
77 use timeCoord_mod
78 use gridStateVector_mod
79 use gridVariableTransforms_mod
80 use bMatrixHI_mod
81 use bMatrixEnsemble_mod
82 use randomNumber_mod
83 use advection_mod
84 use ensembleStateVector_mod
85 use localization_mod
86 use lamBmatrixHI_mod
87 implicit none
88
89 type(struct_vco), pointer :: vco_anl => null()
90 type(struct_hco), pointer :: hco_anl => null()
91 type(struct_hco), pointer :: hco_core => null()
92
93 type(struct_gsv) :: statevector_x
94 type(struct_gsv) :: statevector_y
95 type(struct_gsv) :: statevector_Ly
96 type(struct_gsv) :: statevector_LTx
97
98 real(8) :: innerProduct1_local, innerProduct2_local, innerProduct1_global, innerProduct2_global
99
100 real(8), allocatable :: controlVector1(:)
101 real(8), allocatable :: controlVector2(:)
102
103 integer :: get_max_rss, ierr, cvDim, nulnam, fnom, fclos
104 character(len=20) :: test ! adjoint test type ('Bhi','Bens','advEns','advGSV','loc')
105 NAMELIST /NAMADT/test
106
107 call ver_printNameAndVersion('adjointTest','Various tests of adjoint codes')
108
109 !
110 !- 1. Settings and module initializations
111 !
112 write(*,*)
113 write(*,*) '> midas-adjointTest: setup - START'
114
115 !- 1.1 mpi
116 call mmpi_initialize
117
118 !- 1.2 timings
119 call tmg_init(mmpi_myid, 'TMG_INFO')
120 call utl_tmg_start(0,'Main')
121
122 !- 1.3 RAM disk usage
123 call ram_setup
124
125 !- 1.4 Temporal grid and set dateStamp from env variable
126 call tim_setup()
127 if (tim_getDateStamp()==0) then
128 call utl_abort('midas-adjointTest: dateStamp was not set')
129 end if
130
131 !- 1.6 Constants
132 if ( mmpi_myid == 0 ) then
133 call mpc_printConstants(6)
134 call pre_printPrecisions
135 end if
136
137 !- 1.8 Variables of the model states
138 call gsv_setup
139
140 !- 1.9 Set the horizontal domain
141 call hco_SetupFromFile(hco_anl, './analysisgrid', 'ANALYSIS') ! IN
142 if ( hco_anl % global ) then
143 hco_core => hco_anl
144 else
145 !- Iniatilized the core (Non-Exteded) analysis grid
146 call hco_SetupFromFile( hco_core, './analysisgrid', 'COREGRID') ! IN
147 end if
148
149 !- 1.10 Initialize the vertical coordinate
150 call vco_SetupFromFile( vco_anl, & ! OUT
151 './analysisgrid') ! IN
152
153 write(*,*)
154 write(*,*) '> midas-adjointTest: setup - END'
155 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
156
157 !- 1.11 Variable transforms
158 call gvt_Setup(hco_anl, hco_core, vco_anl)
159
160 !- 1.12 Test selection
161 test = 'Bhi' ! default test
162
163 if ( .not. utl_isNamelistPresent('NAMADT','./flnml') ) then
164 if (mmpi_myid == 0) then
165 write(*,*) 'midas-adjointTest: namadt is missing in the namelist. '&
166 //'The default values will be taken.'
167 end if
168 else
169 nulnam=0
170 ierr=fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
171 read(nulnam, nml=namadt, iostat=ierr)
172 if(ierr.ne.0) call utl_abort('midas-adjointTest: Error reading namelist')
173 if( mmpi_myid == 0 ) write(*,nml=namadt)
174 ierr=fclos(nulnam)
175 end if
176
177 !
178 !- 2. The tests
179 !
180
181 write(*,*)
182 write(*,*) '> midas-adjointTest: '//test
183 if ( test == 'Bhi') then
184 !- 2.1 Bhi
185 call check_bhi
186 else if ( test == 'Bens') then
187 !- 2.2 Bens
188 call check_bens
189 else if ( test == 'advEns') then
190 !- 2.3 AdvectionENS
191 call check_advectionENS
192 else if ( test == 'advGSV') then
193 !- 2.4 AdvectionGSV
194 call check_advectionGSV
195 else if ( test == 'loc') then
196 !- 2.5 Localization
197 call check_loc
198 !else if ( test == 'addMem') then
199 ! !- 2.6 Localization
200 ! call check_addmem
201 else
202 call utl_abort('midas-adjointTest: inexistant test label ('//test//')')
203 end if
204
205 !
206 !- 3. Ending
207 !
208 write(*,*)
209 write(*,*) '> midas-adjointTest: Ending'
210 call utl_tmg_stop(0)
211 call tmg_terminate(mmpi_myid, 'TMG_INFO')
212
213 call rpn_comm_finalize(ierr)
214
215contains
216
217 !--------------------------------------------------------------------------
218 !- check Bhi
219 !--------------------------------------------------------------------------
220 subroutine check_bhi()
221 implicit none
222
223 ! Locals:
224 integer :: seed, kIndex, stepIndex, latIndex, lonIndex, cvIndex
225 real(8), pointer :: field4d_r8(:,:,:,:), field3d_Ly_r8(:,:,:), field3d_x_r8(:,:,:)
226
227 call gvt_setupRefFromTrialFiles('HU')
228 if (hco_anl%global) then
229 call bhi_Setup( hco_anl, vco_anl, & ! IN
230 cvdim ) ! OUT
231 else
232 call lbhi_Setup( hco_anl, hco_core, vco_anl, & ! IN
233 cvdim ) ! OUT
234 end if
235
236 call gsv_allocate(statevector_x , tim_nstepobsinc, hco_anl, vco_anl, &
237 mpi_local_opt=.true., &
238 allocHeight_opt=.false., allocPressure_opt=.false.)
239
240 call gsv_allocate(statevector_Ly , tim_nstepobsinc, hco_anl, vco_anl, &
241 mpi_local_opt=.true., &
242 allocHeight_opt=.false., allocPressure_opt=.false.)
243
244 allocate ( controlVector1(cvDim) )
245
246 ! x
247 !statevector_x%gd3d_r8(:,:,:) = 13.3d0
248 call gsv_getField(statevector_x,field4d_r8)
249 seed=1
250 call rng_setup(abs(seed+mmpi_myid))
251 do kIndex = statevector_x%mykBeg, statevector_x%mykEnd
252 do stepIndex = 1, statevector_x%numStep
253 do latIndex = statevector_x%myLatBeg, statevector_x%myLatEnd
254 do lonIndex = statevector_x%myLonBeg, statevector_x%myLonEnd
255 field4d_r8(lonIndex,latIndex,kIndex,stepIndex) = rng_gaussian()
256 end do
257 end do
258 end do
259 end do
260
261 ! y
262 !controlVector1(:) = 2.4d0
263 do cvIndex = 1, cvDim
264 controlVector1(cvIndex) = rng_gaussian()
265 end do
266
267 ! y
268 !controlVector1 = 0.d0
269 !if (hco_anl%global) then
270 ! call bhi_bSqrtAd( statevector_x, & ! IN
271 ! controlVector1) ! OUT
272 !else
273 ! call lbhi_bSqrtAdj( statevector_x, & ! IN
274 ! controlVector1) ! OUT
275 !end if
276
277 ! Ly
278 if (hco_anl%global) then
279 call bhi_bSqrt( controlVector1, & ! IN
280 statevector_Ly ) ! OUT
281 else
282 call lbhi_bSqrt( controlVector1, & ! IN
283 statevector_Ly ) ! OUT
284 end if
285
286
287 ! <x ,L(y)>
288 innerProduct1_local = 0.d0
289 call gsv_getField(statevector_Ly,field3d_Ly_r8)
290 call gsv_getField(statevector_x, field3d_x_r8 )
291 call euclid(innerProduct1_local, &
292 field3d_x_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:), &
293 field3d_Ly_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:), &
294 statevector_Ly%myLonBeg, statevector_Ly%myLonEnd, statevector_Ly%myLatBeg, statevector_Ly%myLatEnd, statevector_Ly%nk, 1)
295 write(*,*) "<x ,L(y)> local = ",innerProduct1_local
296 call rpn_comm_allreduce(innerProduct1_local,innerProduct1_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
297 write(*,*) "<x ,L(y)> global= ",innerProduct1_global
298
299 ! L_T(x)
300 allocate ( controlVector2(cvDim) )
301 controlVector2 = 0.d0
302 if (hco_anl%global) then
303 call bhi_bSqrtAd( statevector_x, & ! IN
304 controlVector2 ) ! OUT
305 else
306 call lbhi_bSqrtAdj( statevector_x, & ! IN
307 controlVector2 ) ! OUT
308 end if
309
310 ! <L_T(x),y>
311 innerProduct2_local = 0.d0
312 call euclid(innerProduct2_local, controlVector2, controlVector1, 1, cvDim, 1, 1, 1, 1)
313 print*,"<Lt(x) ,y > local = ",innerProduct2_local
314 call rpn_comm_allreduce(innerProduct2_local,innerProduct2_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
315 write(*,*) "<Lt(x) ,y > global= ",innerProduct2_global
316
317 ! Results
318 call checkAndOutputInnerProd
319
320 deallocate(controlVector2)
321 deallocate(controlVector1)
322
323 call gsv_deallocate(statevector_x)
324 call gsv_deallocate(statevector_Ly)
325
326 end subroutine check_bhi
327
328 !--------------------------------------------------------------------------
329 !- check Bens
330 !--------------------------------------------------------------------------
331 subroutine check_bens()
332 implicit none
333
334 ! Locals:
335 integer :: seed, kIndex, stepIndex, latIndex, lonIndex
336 integer, allocatable :: cvDimPerInstance(:)
337 real(8), pointer :: field4d_Ly_r8(:,:,:,:), field4d_x_r8(:,:,:,:)
338
339 call gvt_setupRefFromTrialFiles('HU')
340 call ben_Setup( hco_anl, hco_core, vco_anl, & ! IN
341 cvdimPerInstance ) ! OUT
342
343 cvDim = cvdimPerInstance(1)
344 call gsv_allocate(statevector_x , tim_nstepobsinc, hco_anl, vco_anl, &
345 mpi_local_opt=.true., &
346 allocHeight_opt=.false., allocPressure_opt=.false.)
347 call gsv_allocate(statevector_Ly , tim_nstepobsinc, hco_anl, vco_anl, &
348 mpi_local_opt=.true., &
349 allocHeight_opt=.false., allocPressure_opt=.false.)
350
351 allocate ( controlVector1(cvdim) )
352
353 ! x
354 seed=1
355 call rng_setup(abs(seed+mmpi_myid))
356 call gsv_getField(statevector_x,field4d_x_r8)
357 do kIndex = statevector_x%mykBeg, statevector_x%mykEnd
358 do stepIndex = 1, statevector_x%numStep
359 do latIndex = statevector_x%myLatBeg, statevector_x%myLatEnd
360 do lonIndex = statevector_x%myLonBeg, statevector_x%myLonEnd
361 field4d_x_r8(lonIndex,latIndex,kIndex,stepIndex) = rng_gaussian()
362 end do
363 end do
364 end do
365 end do
366
367 ! y
368 controlVector1 = 0.d0
369 call ben_bSqrtAd( 1, statevector_x, & ! IN
370 controlVector1) ! OUT
371
372 ! Ly
373 call ben_bSqrt( 1, controlVector1, & ! IN
374 statevector_Ly ) ! OUT
375
376 ! <x ,L(y)>
377 innerProduct1_local = 0.d0
378 call gsv_getField(statevector_Ly,field4d_Ly_r8)
379 call gsv_getField(statevector_x, field4d_x_r8 )
380 call euclid(innerProduct1_local, &
381 field4d_x_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:,:), &
382 field4d_Ly_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:,:), &
383 statevector_Ly%myLonBeg, statevector_Ly%myLonEnd, statevector_Ly%myLatBeg, statevector_Ly%myLatEnd, statevector_Ly%nk, statevector_Ly%numStep)
384 write(*,*) "<x ,L(y)> local = ",innerProduct1_local
385 call rpn_comm_allreduce(innerProduct1_local,innerProduct1_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
386 write(*,*) "<x ,L(y)> global= ",innerProduct1_global
387
388 ! L_T(x)
389 allocate ( controlVector2(cvDim) )
390 controlVector2 = 0.d0
391 call ben_bSqrtAd( 1, statevector_x, & ! IN
392 controlVector2 ) ! OUT
393
394 ! <L_T(x),y>
395 innerProduct2_local = 0.d0
396 call euclid(innerProduct2_local, controlVector2, controlVector1, 1, cvDim, 1, 1, 1, 1)
397 print*,"<Lt(x) ,y > local = ",innerProduct2_local
398 call rpn_comm_allreduce(innerProduct2_local,innerProduct2_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
399 write(*,*) "<Lt(x) ,y > global= ",innerProduct2_global
400
401 ! Results
402 call checkAndOutputInnerProd
403
404 deallocate(controlVector2)
405 deallocate(controlVector1)
406
407 call gsv_deallocate(statevector_x)
408 call gsv_deallocate(statevector_Ly)
409
410 end subroutine check_bens
411
412 !--------------------------------------------------------------------------
413 !- check Loc
414 !--------------------------------------------------------------------------
415 subroutine check_loc()
416 implicit none
417
418 ! Locals:
419 integer :: seed, kIndex, stepIndex, latIndex, lonIndex, dateStamp
420 integer :: numStepAmplitude, amp3dStepIndex, memberIndex
421 type(struct_ens) :: ensAmplitude_x
422 type(struct_ens) :: ensAmplitude_Ly
423 type(struct_loc), pointer :: loc => null()
424 real(8), pointer :: ens_oneLev(:,:,:,:)
425 real(8), pointer :: field4d_Ly_r8(:,:,:,:), field4d_x_r8(:,:,:,:)
426 character(len=4), parameter :: varNameALFAatm(1) = (/ 'ALFA' /)
427 character(len=4), parameter :: varNameALFAsfc(1) = (/ 'ALFS' /)
428 character(len=4) :: varNameALFA(1)
429 integer,allocatable :: dateStampList(:)
430 integer, allocatable :: cvDimPerInstance(:)
431
432 allocate(dateStampList(tim_nstepobsinc))
433 call tim_getstamplist(dateStampList,tim_nstepobsinc,tim_getDatestamp())
434
435 dateStamp = tim_getDatestamp()
436 write(*,*) 'check_loc: tim_getDatestamp = ', dateStamp
437 write(*,*) 'check_loc: dateStampList = ', dateStampList(:)
438
439 call ben_Setup( hco_anl, hco_core, vco_anl, & ! IN
440 cvDimPerInstance ) ! OUT
441
442 cvDim = cvDimPerInstance(1)
443
444 loc => ben_getLoc(1)
445 if ( cvDim /= loc%cvDim ) then
446 call utl_abort('check_loc: cvDim /= loc%cvDim')
447 end if
448 numStepAmplitude = ben_getNumStepAmplitudeAssimWindow()
449 if ( numStepAmplitude /= 1 ) then
450 call utl_abort('check_loc: not yet adapted for localization advection')
451 end if
452 amp3dStepIndex = ben_getAmp3dStepIndexAssimWindow()
453
454 if (loc%vco%Vcode == 5002 .or. loc%vco%Vcode == 5005) then
455 varNameALFA(:) = varNameALFAatm(:)
456 else ! vco_anl%Vcode == -1
457 varNameALFA(:) = varNameALFAsfc(:)
458 end if
459
460 call gsv_allocate(statevector_x , numStepAmplitude, loc%hco, loc%vco, &
461 mpi_local_opt=.true., varNames_opt=varNameALFA, dataKind_opt=8)
462 call gsv_allocate(statevector_Ly , numStepAmplitude, loc%hco, loc%vco, &
463 mpi_local_opt=.true., varNames_opt=varNameALFA, dataKind_opt=8)
464
465 call ens_allocate(ensAmplitude_x, loc%nEnsOverDimension, numStepAmplitude, loc%hco, loc%vco, &
466 datestampList=dateStampList, varNames_opt=varNameALFA, dataKind_opt=8)
467 call ens_allocate(ensAmplitude_Ly, loc%nEnsOverDimension, numStepAmplitude, loc%hco, loc%vco, &
468 datestampList=dateStampList, varNames_opt=varNameALFA, dataKind_opt=8)
469
470 allocate ( controlVector1(cvDim) )
471
472 ! x
473 seed=1
474 call rng_setup(abs(seed+mmpi_myid))
475 do kIndex = 1, ens_getNumK(ensAmplitude_x)
476 ens_oneLev => ens_getOneLev_r8(ensAmplitude_x,kIndex)
477 do memberIndex = 1, loc%nEnsOverDimension
478 do stepIndex = 1,numStepAmplitude
479 do latIndex = statevector_x%myLatBeg, statevector_x%myLatEnd
480 do lonIndex = statevector_x%myLonBeg, statevector_x%myLonEnd
481 ens_oneLev(memberIndex,stepIndex,lonIndex,latIndex) = rng_gaussian()
482 end do
483 end do
484 end do
485 end do
486 end do
487
488 ! y
489 controlVector1 = 0.d0
490 call loc_LsqrtAd(loc, & ! IN
491 ensAmplitude_x, & ! IN
492 controlVector1, & ! OUT
493 amp3dStepIndex) ! IN
494
495 ! Ly
496 call loc_Lsqrt (loc, & ! IN
497 controlVector1, & ! IN
498 ensAmplitude_Ly, & ! OUT
499 amp3dStepIndex) ! IN
500
501 ! <x ,L(y)>
502 innerProduct1_local = 0.d0
503 call gsv_getField(statevector_Ly,field4d_Ly_r8)
504 call gsv_getField(statevector_x, field4d_x_r8 )
505 do memberIndex = 1, loc%nEnsOverDimension
506 call ens_copyMember(ensAmplitude_x , statevector_x , memberIndex)
507 call ens_copyMember(ensAmplitude_Ly, statevector_Ly, memberIndex)
508 call euclid(innerProduct1_local, &
509 field4d_x_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:,:), &
510 field4d_Ly_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:,:), &
511 statevector_Ly%myLonBeg, statevector_Ly%myLonEnd, statevector_Ly%myLatBeg, statevector_Ly%myLatEnd, statevector_Ly%nk, statevector_Ly%numStep)
512 end do
513 write(*,*) "<x ,L(y)> local = ",innerProduct1_local
514 call rpn_comm_allreduce(innerProduct1_local,innerProduct1_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
515 write(*,*) "<x ,L(y)> global= ",innerProduct1_global
516
517 ! L_T(x)
518 allocate ( controlVector2(cvDim) )
519 controlVector2 = 0.d0
520 call loc_LsqrtAd(loc, & ! IN
521 ensAmplitude_x, & ! IN
522 controlVector2, & ! OUT
523 amp3dStepIndex) ! IN
524
525 ! <L_T(x),y>
526 innerProduct2_local = 0.d0
527 call euclid(innerProduct2_local, controlVector2, controlVector1, 1, cvDim, 1, 1, 1, 1)
528 print*,"<Lt(x) ,y > local = ",innerProduct2_local
529 call rpn_comm_allreduce(innerProduct2_local,innerProduct2_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
530 write(*,*) "<Lt(x) ,y > global= ",innerProduct2_global
531
532 ! Results
533 call checkAndOutputInnerProd
534
535 deallocate(controlVector2)
536 deallocate(controlVector1)
537
538 call gsv_deallocate(statevector_x)
539 call gsv_deallocate(statevector_Ly)
540
541 call ens_deallocate(ensAmplitude_x)
542 call ens_deallocate(ensAmplitude_Ly)
543
544 end subroutine check_loc
545
546!!$ !--------------------------------------------------------------------------
547!!$ !- check addMem
548!!$ !--------------------------------------------------------------------------
549!!$ subroutine check_addmem()
550!!$ implicit none
551!!$
552!!$ integer :: seed, kIndex, stepIndex, latIndex, lonIndex
553!!$ integer :: numStepAmplitude, amp3dStepIndex, memberIndex
554!!$
555!!$ type(struct_ens) :: ensAmplitude_LTx
556!!$ type(struct_ens) :: ensAmplitude_y
557!!$
558!!$ type(struct_loc), pointer :: loc => null()
559!!$
560!!$ real(8), pointer :: ens_oneLev(:,:,:,:)
561!!$
562!!$ character(len=4), parameter :: varNameALFAatm(1) = (/ 'ALFA' /)
563!!$ character(len=4), parameter :: varNameALFAsfc(1) = (/ 'ALFS' /)
564!!$ character(len=4) :: varNameALFA(1)
565!!$
566!!$ integer,allocatable :: dateStampList(:)
567!!$
568!!$ allocate(dateStampList(tim_nstepobsinc))
569!!$ call tim_getstamplist(dateStampList,tim_nstepobsinc,tim_getDatestamp())
570!!$
571!!$ write(*,*) 'JFC tim_getDatestamp = ', tim_getDatestamp()
572!!$ write(*,*) 'JFC dateStampList = ', dateStampList(:)
573!!$
574!!$ call ben_Setup( hco_anl, hco_core, vco_anl, & ! IN
575!!$ cvdim ) ! OUT
576!!$
577!!$ write(*,*) 'JFC ben_Setup done '
578!!$
579!!$
580!!$ loc => ben_getLoc(1)
581!!$ if ( cvDim /= loc%cvDim ) then
582!!$ call utl_abort('check_loc: cvDim /= loc%cvDim')
583!!$ end if
584!!$ numStepAmplitude = ben_getNumStepAmplitudeAssimWindow()
585!!$ if ( numStepAmplitude /= 1 ) then
586!!$ call utl_abort('check_loc: not yet adapted for localization advection')
587!!$ end if
588!!$ amp3dStepIndex = ben_getAmp3dStepIndexAssimWindow()
589!!$
590!!$ if (loc%vco%Vcode == 5002 .or. loc%vco%Vcode == 5005) then
591!!$ varNameALFA(:) = varNameALFAatm(:)
592!!$ else ! vco_anl%Vcode == -1
593!!$ varNameALFA(:) = varNameALFAsfc(:)
594!!$ end if
595!!$
596!!$ call gsv_allocate(statevector_x , tim_nstepobsinc, hco_anl, vco_anl, &
597!!$ mpi_local_opt=.true., &
598!!$ allocHeight_opt=.false., allocPressure_opt=.false.)
599!!$ call gsv_allocate(statevector_Ly , tim_nstepobsinc, hco_anl, vco_anl, &
600!!$ mpi_local_opt=.true., &
601!!$ allocHeight_opt=.false., allocPressure_opt=.false.)
602!!$
603!!$ call gsv_allocate(statevector_LTx , numStepAmplitude, loc%hco, loc%vco, &
604!!$ mpi_local_opt=.true., varNames_opt=varNameALFA, dataKind_opt=8)
605!!$ call gsv_allocate(statevector_y , numStepAmplitude, loc%hco, loc%vco, &
606!!$ mpi_local_opt=.true., varNames_opt=varNameALFA, dataKind_opt=8)
607!!$
608!!$ call ens_allocate(ensAmplitude_LTx, loc%nEnsOverDimension, numStepAmplitude, loc%hco, loc%vco, &
609!!$ datestampList=dateStampList, varNames_opt=varNameALFA, dataKind_opt=8)
610!!$ call ens_allocate(ensAmplitude_y, loc%nEnsOverDimension, numStepAmplitude, loc%hco, loc%vco, &
611!!$ datestampList=dateStampList, varNames_opt=varNameALFA, dataKind_opt=8)
612!!$
613!!$ ! x
614!!$ seed=1
615!!$ call rng_setup(abs(seed+mmpi_myid))
616!!$ do kIndex = statevector_x%mykBeg, statevector_x%mykEnd
617!!$ do stepIndex = 1, statevector_x%numStep
618!!$ do latIndex = statevector_x%myLatBeg, statevector_x%myLatEnd
619!!$ do lonIndex = statevector_x%myLonBeg, statevector_x%myLonEnd
620!!$ statevector_x%gd_r8(lonIndex,latIndex,kIndex,stepIndex) = rng_gaussian()
621!!$ end do
622!!$ end do
623!!$ end do
624!!$ end do
625!!$
626!!$ ! y
627!!$ do kIndex = 1, ens_getNumK(ensAmplitude_y)
628!!$ ens_oneLev => ens_getOneLev_r8(ensAmplitude_y,kIndex)
629!!$ do memberIndex = 1, loc%nEnsOverDimension
630!!$ do stepIndex = 1,tim_nstepobsinc
631!!$ do latIndex = statevector_x%myLatBeg, statevector_x%myLatEnd
632!!$ do lonIndex = statevector_x%myLonBeg, statevector_x%myLonEnd
633!!$ ens_oneLev(memberIndex,stepIndex,lonIndex,latIndex) = rng_gaussian()
634!!$ end do
635!!$ end do
636!!$ end do
637!!$ end do
638!!$ end do
639!!$
640!!$ ! Ly
641!!$ call addEnsMember( ensAmplitude_y, & ! IN
642!!$ statevector_Ly, & ! OUT
643!!$ 1, .false. ) ! IN
644!!$
645!!$ ! <x ,L(y)>
646!!$ innerProduct1_local = 0.d0
647!!$ call euclid(innerProduct1_local, &
648!!$ statevector_x %gd_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:,:), &
649!!$ statevector_Ly%gd_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:,:), &
650!!$ statevector_Ly%myLonBeg, statevector_Ly%myLonEnd, statevector_Ly%myLatBeg, statevector_Ly%myLatEnd, statevector_Ly%nk, statevector_Ly%numStep)
651!!$ write(*,*) "<x ,L(y)> local = ",innerProduct1_local
652!!$ call rpn_comm_allreduce(innerProduct1_local,innerProduct1_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
653!!$ write(*,*) "<x ,L(y)> global= ",innerProduct1_global
654!!$
655!!$ ! L_T(x)
656!!$ call addEnsMemberad( statevector_x, & ! IN
657!!$ ensAmplitude_LTx, & ! OUT
658!!$ 1, .false. ) ! IN
659!!$
660!!$ ! <L_T(x),y>
661!!$ innerProduct2_local = 0.d0
662!!$ do memberIndex = 1, loc%nEnsOverDimension
663!!$ call ens_copyMember(ensAmplitude_LTx, statevector_LTx, memberIndex)
664!!$ call ens_copyMember(ensAmplitude_y , statevector_y , memberIndex)
665!!$ call euclid(innerProduct2_local, &
666!!$ statevector_LTx %gd_r8(statevector_y%myLonBeg:statevector_y%myLonEnd,statevector_y%myLatBeg:statevector_y%myLatEnd,:,:), &
667!!$ statevector_y%gd_r8(statevector_y%myLonBeg:statevector_y%myLonEnd,statevector_y%myLatBeg:statevector_y%myLatEnd,:,:), &
668!!$ statevector_y%myLonBeg, statevector_y%myLonEnd, statevector_y%myLatBeg, statevector_y%myLatEnd, statevector_y%nk, statevector_y%numStep)
669!!$ end do
670!!$ print*,"<Lt(x) ,y > local = ",innerProduct2_local
671!!$ call rpn_comm_allreduce(innerProduct2_local,innerProduct2_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
672!!$ write(*,*) "<Lt(x) ,y > global= ",innerProduct2_global
673!!$
674!!$ ! Results
675!!$ call checkAndOutputInnerProd
676!!$
677!!$ call gsv_deallocate(statevector_LTx)
678!!$ call gsv_deallocate(statevector_y)
679!!$ call gsv_deallocate(statevector_x)
680!!$ call gsv_deallocate(statevector_Ly)
681!!$
682!!$ call ens_deallocate(ensAmplitude_LTx)
683!!$ call ens_deallocate(ensAmplitude_y)
684!!$
685!!$ end subroutine check_addmem
686
687 !--------------------------------------------------------------------------
688 !- check AdvectionENS
689 !--------------------------------------------------------------------------
690 subroutine check_advectionENS()
691 implicit none
692
693 ! Locals:
694 integer :: seed, kIndex, stepIndex, latIndex, lonIndex, dateStamp
695 type(struct_adv) :: adv_analInc
696 type(struct_ens) :: ens_x, ens_Ly, ens_y, ens_LTx
697 character(len=32) :: directionAnlInc
698 character(len=4), parameter :: varNameALFA(1) = (/ 'ALFA' /)
699 real(8), pointer :: ens_oneLev(:,:,:,:)
700 real(8), pointer :: field4d_Ly_r8(:,:,:,:), field4d_x_r8(:,:,:,:)
701 real(8), pointer :: field4d_LTx_r8(:,:,:,:), field4d_y_r8(:,:,:,:)
702 real(8) :: delT_hour
703 real(8), allocatable :: advectFactor(:)
704 integer :: numStepAdvect, numStepReferenceFlow, nEns, memberIndex
705 integer,allocatable :: dateStampList(:)
706
707 allocate(dateStampList(tim_nstepobsinc))
708 call tim_getstamplist(dateStampList,tim_nstepobsinc,tim_getDatestamp())
709
710 dateStamp = tim_getDatestamp()
711 write(*,*) 'check_advectionENS: tim_getDatestamp = ', dateStamp
712 write(*,*) 'check_advectionENS: dateStampList = ', dateStampList(:)
713
714 directionAnlInc = 'fromFirstTimeIndex' !'towardFirstTimeIndexInverse'
715 delT_hour = 1.0d0 !tim_dstepobsinc
716 numStepAdvect = tim_nstepobsinc
717 numStepReferenceFlow = 7
718 allocate(advectFactor(vco_anl%nLev_M))
719 advectFactor(:) = 0.75D0
720 nEns = 1
721
722 call adv_setup( adv_analInc, & ! OUT
723 directionAnlInc, hco_anl, vco_anl, & ! IN
724 numStepAdvect, dateStampList, & ! IN
725 numStepReferenceFlow, delT_hour, advectFactor, & ! IN
726 'MMLevsOnly', steeringFlowFilename_opt='ensemble/forecast_for_advection' ) ! IN
727
728 deallocate(advectFactor)
729
730 call ens_allocate(ens_x, nEns, numStepAdvect, hco_anl, vco_anl, dateStampList, &
731 hco_core_opt=hco_core, varNames_opt=varNameALFA, dataKind_opt=8)
732 call ens_allocate(ens_Ly, nEns, numStepAdvect, hco_anl, vco_anl, dateStampList, &
733 hco_core_opt=hco_core, varNames_opt=varNameALFA, dataKind_opt=8)
734 call ens_allocate(ens_y, nEns, numStepAdvect, hco_anl, vco_anl, dateStampList, &
735 hco_core_opt=hco_core, varNames_opt=varNameALFA, dataKind_opt=8)
736 call ens_allocate(ens_LTx, nEns, numStepAdvect, hco_anl, vco_anl, dateStampList, &
737 hco_core_opt=hco_core, varNames_opt=varNameALFA, dataKind_opt=8)
738
739 call gsv_allocate(statevector_x , tim_nstepobsinc, hco_anl, vco_anl, &
740 mpi_local_opt=.true., &
741 allocHeight_opt=.false., allocPressure_opt=.false.)
742 call gsv_allocate(statevector_Ly , tim_nstepobsinc, hco_anl, vco_anl, &
743 mpi_local_opt=.true., &
744 allocHeight_opt=.false., allocPressure_opt=.false.)
745 call gsv_allocate(statevector_y , tim_nstepobsinc, hco_anl, vco_anl, &
746 mpi_local_opt=.true., &
747 allocHeight_opt=.false., allocPressure_opt=.false.)
748 call gsv_allocate(statevector_LTx , tim_nstepobsinc, hco_anl, vco_anl, &
749 mpi_local_opt=.true., &
750 allocHeight_opt=.false., allocPressure_opt=.false.)
751
752 ! x
753 seed=1
754 call rng_setup(abs(seed+mmpi_myid))
755 do kIndex = 1, ens_getNumK(ens_x)
756 ens_oneLev => ens_getOneLev_r8(ens_x,kIndex)
757 do memberIndex = 1, nEns
758 do stepIndex = 1,tim_nstepobsinc
759 do latIndex = statevector_x%myLatBeg, statevector_x%myLatEnd
760 do lonIndex = statevector_x%myLonBeg, statevector_x%myLonEnd
761 ens_oneLev(memberIndex,stepIndex,lonIndex,latIndex) = rng_gaussian()
762 end do
763 end do
764 end do
765 end do
766 end do
767
768 ! y
769 call ens_copy(ens_x, & ! IN
770 ens_y) ! OUT
771 call adv_ensemble_ad( ens_y, & ! INOUT
772 adv_analInc, nEns ) ! IN
773 ! Ly
774 call ens_copy(ens_y, ens_Ly)
775 call adv_ensemble_tl( ens_Ly, & ! INOUT
776 adv_analInc, nEns ) ! IN
777
778 ! <x ,L(y)>
779 innerProduct1_local = 0.d0
780 call gsv_getField(statevector_Ly,field4d_Ly_r8)
781 call gsv_getField(statevector_x, field4d_x_r8 )
782 do memberIndex = 1, nEns
783 call ens_copyMember(ens_x , statevector_x , memberIndex)
784 call ens_copyMember(ens_Ly, statevector_Ly, memberIndex)
785 call euclid(innerProduct1_local, &
786 field4d_x_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:,:), &
787 field4d_Ly_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:,:), &
788 statevector_Ly%myLonBeg, statevector_Ly%myLonEnd, statevector_Ly%myLatBeg, statevector_Ly%myLatEnd, statevector_Ly%nk, statevector_Ly%numStep)
789 end do
790 write(*,*) "<x ,L(y)> local = ",innerProduct1_local
791 call rpn_comm_allreduce(innerProduct1_local,innerProduct1_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
792 write(*,*) "<x ,L(y)> global= ",innerProduct1_global
793
794 ! L_T(x)
795 call ens_copy(ens_x, & ! IN
796 ens_LTx) ! OUT
797 call adv_ensemble_ad( ens_LTx, & ! INOUT
798 adv_analInc, nEns ) ! IN
799
800 ! <L_T(x),y>
801 innerProduct2_local = 0.d0
802 call gsv_getField(statevector_LTx,field4d_LTx_r8)
803 call gsv_getField(statevector_y, field4d_y_r8 )
804 do memberIndex = 1, nEns
805 call ens_copyMember(ens_LTx, statevector_LTx, memberIndex)
806 call ens_copyMember(ens_y , statevector_y , memberIndex)
807 call euclid(innerProduct2_local, &
808 field4d_LTx_r8(statevector_y%myLonBeg:statevector_y%myLonEnd,statevector_y%myLatBeg:statevector_y%myLatEnd,:,:), &
809 field4d_y_r8(statevector_y%myLonBeg:statevector_y%myLonEnd,statevector_y%myLatBeg:statevector_y%myLatEnd,:,:), &
810 statevector_y%myLonBeg, statevector_y%myLonEnd, statevector_y%myLatBeg, statevector_y%myLatEnd, statevector_y%nk, statevector_y%numStep)
811 end do
812 print*,"<Lt(x) ,y > local = ",innerProduct2_local
813 call rpn_comm_allreduce(innerProduct2_local,innerProduct2_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
814 write(*,*) "<Lt(x) ,y > global= ",innerProduct2_global
815
816 ! Results
817 call checkAndOutputInnerProd
818
819 call gsv_deallocate(statevector_x)
820 call gsv_deallocate(statevector_Ly)
821 call gsv_deallocate(statevector_LTx)
822 call gsv_deallocate(statevector_y)
823
824 call ens_deallocate(ens_x)
825 call ens_deallocate(ens_Ly)
826 call ens_deallocate(ens_LTx)
827 call ens_deallocate(ens_y)
828
829 end subroutine check_advectionENS
830
831 !--------------------------------------------------------------------------
832 !- check AdvectionGSV
833 !--------------------------------------------------------------------------
834 subroutine check_advectionGSV()
835 implicit none
836
837 ! Locals:
838 integer :: seed, kIndex, stepIndex, latIndex, lonIndex, dateStamp
839 type(struct_adv) :: adv_analInc
840 character(len=32) :: directionAnlInc
841 real(8) :: delT_hour
842 real(8), allocatable :: advectFactor(:)
843 real(8), pointer :: field4d_x_r8(:,:,:,:), field4d_y_r8(:,:,:,:)
844 real(8), pointer :: field4d_LTx_r8(:,:,:,:), field4d_Ly_r8(:,:,:,:)
845 integer :: numStepAdvect, numStepReferenceFlow
846 integer,allocatable :: dateStampList(:)
847
848 allocate(dateStampList(tim_nstepobsinc))
849 call tim_getstamplist(dateStampList,tim_nstepobsinc,tim_getDatestamp())
850
851 dateStamp = tim_getDatestamp()
852 write(*,*) 'check_advectionGSV: tim_getDatestamp = ', dateStamp
853 write(*,*) 'check_advectionGSV: dateStampList = ', dateStampList(:)
854
855 directionAnlInc = 'fromFirstTimeIndex' !'towardFirstTimeIndexInverse'
856 delT_hour = 1.0d0 !tim_dstepobsinc
857 numStepAdvect = tim_nstepobsinc
858 numStepReferenceFlow = 7
859 advectFactor = 0.75D0
860
861 call adv_setup( adv_analInc, & ! OUT
862 directionAnlInc, hco_anl, vco_anl, & ! IN
863 numStepAdvect, dateStampList, & ! IN
864 numStepReferenceFlow, delT_hour, advectFactor, & ! IN
865 'allLevs', steeringFlowFilename_opt='ensemble/forecast_for_advection' ) ! IN
866
867 call gsv_allocate(statevector_x , tim_nstepobsinc, hco_anl, vco_anl, &
868 mpi_local_opt=.true., &
869 allocHeight_opt=.false., allocPressure_opt=.false.)
870 call gsv_allocate(statevector_Ly , tim_nstepobsinc, hco_anl, vco_anl, &
871 mpi_local_opt=.true., &
872 allocHeight_opt=.false., allocPressure_opt=.false.)
873 call gsv_allocate(statevector_y , tim_nstepobsinc, hco_anl, vco_anl, &
874 mpi_local_opt=.true., &
875 allocHeight_opt=.false., allocPressure_opt=.false.)
876 call gsv_allocate(statevector_LTx , tim_nstepobsinc, hco_anl, vco_anl, &
877 mpi_local_opt=.true., &
878 allocHeight_opt=.false., allocPressure_opt=.false.)
879
880 ! x
881 seed=1
882 call rng_setup(abs(seed+mmpi_myid))
883 call gsv_getField(statevector_x, field4d_x_r8 )
884 do kIndex = statevector_x%mykBeg, statevector_x%mykEnd
885 do stepIndex = 1, statevector_x%numStep
886 do latIndex = statevector_x%myLatBeg, statevector_x%myLatEnd
887 do lonIndex = statevector_x%myLonBeg, statevector_x%myLonEnd
888 field4d_x_r8(lonIndex,latIndex,kIndex,stepIndex) = rng_gaussian()
889 end do
890 end do
891 end do
892 end do
893
894 ! y
895 call gsv_getField(statevector_y, field4d_y_r8 )
896 do kIndex = statevector_y%mykBeg, statevector_y%mykEnd
897 do stepIndex = 1, statevector_y%numStep
898 do latIndex = statevector_y%myLatBeg, statevector_y%myLatEnd
899 do lonIndex = statevector_y%myLonBeg, statevector_y%myLonEnd
900 field4d_y_r8(lonIndex,latIndex,kIndex,stepIndex) = rng_gaussian()
901 end do
902 end do
903 end do
904 end do
905
906 ! Ly
907 call gsv_copy(statevector_y, & ! IN
908 statevector_Ly) ! OUT
909 call adv_statevector_tl( statevector_Ly, & ! INOUT
910 adv_analInc ) ! IN
911
912 ! <x ,L(y)>
913 innerProduct1_local = 0.d0
914 call gsv_getField(statevector_x, field4d_x_r8 )
915 call gsv_getField(statevector_Ly, field4d_Ly_r8 )
916 call euclid(innerProduct1_local, &
917 field4d_x_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:,:), &
918 field4d_Ly_r8(statevector_Ly%myLonBeg:statevector_Ly%myLonEnd,statevector_Ly%myLatBeg:statevector_Ly%myLatEnd,:,:), &
919 statevector_Ly%myLonBeg, statevector_Ly%myLonEnd, statevector_Ly%myLatBeg, statevector_Ly%myLatEnd, statevector_Ly%nk, statevector_Ly%numStep)
920 write(*,*) "<x ,L(y)> local = ",innerProduct1_local
921 call rpn_comm_allreduce(innerProduct1_local,innerProduct1_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
922 write(*,*) "<x ,L(y)> global= ",innerProduct1_global
923
924 ! L_T(x)
925 call gsv_copy(statevector_x, & ! IN
926 statevector_LTx) ! OUT
927 call adv_statevector_ad( statevector_LTx, & ! INOUT
928 adv_analInc ) ! IN
929
930 ! <L_T(x),y>
931 innerProduct2_local = 0.d0
932 call gsv_getField(statevector_LTx, field4d_LTx_r8 )
933 call gsv_getField(statevector_y, field4d_y_r8 )
934 call euclid(innerProduct2_local, &
935 field4d_LTx_r8(statevector_y%myLonBeg:statevector_y%myLonEnd,statevector_y%myLatBeg:statevector_y%myLatEnd,:,:), &
936 field4d_y_r8(statevector_y%myLonBeg:statevector_y%myLonEnd,statevector_y%myLatBeg:statevector_y%myLatEnd,:,:), &
937 statevector_y%myLonBeg, statevector_y%myLonEnd, statevector_y%myLatBeg, statevector_y%myLatEnd, statevector_y%nk, statevector_y%numStep)
938! call euclid(innerProduct2_local, controlVector2, controlVector1, 1, cvDim, 1, 1, 1, 1)
939 print*,"<Lt(x) ,y > local = ",innerProduct2_local
940 call rpn_comm_allreduce(innerProduct2_local,innerProduct2_global,1,"mpi_double_precision","mpi_sum","GRID",ierr)
941 write(*,*) "<Lt(x) ,y > global= ",innerProduct2_global
942
943 ! Results
944 call checkAndOutputInnerProd
945
946 call gsv_deallocate(statevector_x)
947 call gsv_deallocate(statevector_Ly)
948 call gsv_deallocate(statevector_LTx)
949 call gsv_deallocate(statevector_y)
950
951 end subroutine check_advectionGSV
952
953 !--------------------------------------------------------------------------
954 !- Inner product computation
955 !--------------------------------------------------------------------------
956 subroutine euclid (pvalue,px,py,nis,nie,njs,nje,nk,nStep)
957 implicit none
958
959 ! Arguments:
960 integer, intent(in) :: nis
961 integer, intent(in) :: nie
962 integer, intent(in) :: njs
963 integer, intent(in) :: nje
964 integer, intent(in) :: nk
965 integer, intent(in) :: nStep
966 real(8), intent(in) :: px(nis:nie,njs:nje,nk,nStep)
967 real(8), intent(in) :: py(nis:nie,njs:nje,nk,nStep)
968 real(8), intent(inout) :: pvalue
969
970 ! Locals:
971 integer :: i,j,k,s
972
973 do s = 1, nStep
974 do k = 1,nk
975 do j = njs, nje
976 do i = nis, nie
977 pvalue = pvalue + px(i,j,k,s) * py(i,j,k,s)
978 end do
979 end do
980 end do
981 end do
982
983 end subroutine euclid
984
985 !--------------------------------------------------------------------------
986 !- Inner product comparison
987 !--------------------------------------------------------------------------
988 subroutine checkAndOutputInnerProd()
989 implicit none
990
991 ! Locals:
992 integer :: fun
993
994 if ( mmpi_myid == 0 ) then
995 open(newunit=fun, file="innerProd.txt", status="new", action="write")
996 write(fun,'(A20)') test
997 write(fun, '(G23.16)') innerProduct1_global
998 write(fun, '(G23.16)') innerProduct2_global
999 write(fun, '(G23.16)') innerProduct2_global-innerProduct1_global
1000 write(*,*)
1001 if ( innerProduct2_global + innerProduct1_global /= 0.d0 ) then
1002 write(fun, '(G23.16,A1)') 100.d0 * (abs(innerProduct2_global-innerProduct1_global) &
1003 /(0.5d0*(innerProduct2_global+innerProduct1_global)) ), &
1004 '%'
1005 write(*,'(A20,1X,A,1X,G23.16)') test, ": InnerProd Difference(%) = ", &
1006 100.d0 * (abs(innerProduct2_global-innerProduct1_global) / &
1007 (0.5d0*(innerProduct2_global+innerProduct1_global)) )
1008 else
1009 write(fun, '(A23)') 'InnerProduct = 0 ERROR'
1010 write(*,*) 'InnerProduct = 0 !!! Obviously, something went wrong...'
1011 end if
1012 close(fun)
1013 end if
1014
1015 end subroutine checkAndOutputInnerProd
1016
1017end program midas_adjointTest
