FV3 Bundle
type_nicas.F90
Go to the documentation of this file.
1 !----------------------------------------------------------------------
2 ! Module: type_nicas
3 ! Purpose: NICAS data derived type
4 ! Author: Benjamin Menetrier
5 ! Licensing: this code is distributed under the CeCILL-C license
6 ! Copyright © 2015-... UCAR, CERFACS, METEO-FRANCE and IRIT
7 !----------------------------------------------------------------------
8 module type_nicas
9 
10 use netcdf
11 use tools_const, only: rad2deg,reqkm,pi,msvali,msvalr
12 use tools_func, only: cholesky
13 use tools_kinds, only: kind_real
14 use tools_missing, only: msi,msr,isnotmsr
15 use tools_nc, only: ncfloat
16 use type_bpar, only: bpar_type
17 use type_cmat, only: cmat_type
18 use type_com, only: com_type
19 use type_cv, only: cv_type
20 use type_ens, only: ens_type
21 use type_geom, only: geom_type
22 use type_hdiag, only: hdiag_type
23 use type_io, only: io_type
24 use type_linop, only: linop_type
25 use type_nicas_blk, only: nicas_blk_type
26 use type_mpl, only: mpl_type
27 use type_nam, only: nam_type
28 use type_rng, only: rng_type
29 use fckit_mpi_module, only: fckit_mpi_sum
30 
31 implicit none
32 
33 integer,parameter :: ne_rand = 150 ! Ensemble size for randomization
34 integer,parameter :: nfac = 10 ! Number of length-scale factors
35 integer,parameter :: ntest = 100 ! Number of tests
36 logical,parameter :: pos_def_test = .false. ! Positive-definiteness test
37 
38 ! NICAS derived type
39 type nicas_type
40  character(len=1024) :: prefix ! Prefix
41  type(nicas_blk_type),allocatable :: blk(:) ! NICAS data blocks
42  logical :: allocated ! Allocation flag
43 contains
44  procedure :: alloc => nicas_alloc
45  procedure :: dealloc => nicas_dealloc
46  procedure :: read => nicas_read
47  procedure :: write => nicas_write
48  procedure :: write_mpi_summary => nicas_write_mpi_summary
49  procedure :: run_nicas => nicas_run_nicas
50  procedure :: run_nicas_tests => nicas_run_nicas_tests
51  procedure :: alloc_cv => nicas_alloc_cv
52  procedure :: random_cv => nicas_random_cv
53  procedure :: apply => nicas_apply
54  procedure :: apply_from_sqrt => nicas_apply_from_sqrt
55  procedure :: apply_sqrt => nicas_apply_sqrt
56  procedure :: apply_sqrt_ad => nicas_apply_sqrt_ad
57  procedure :: randomize => nicas_randomize
58  procedure :: apply_bens => nicas_apply_bens
59  procedure :: apply_bens_noloc => nicas_apply_bens_noloc
60  procedure :: test_adjoint => nicas_test_adjoint
61  procedure :: test_sqrt => nicas_test_sqrt
62  procedure :: test_dirac => nicas_test_dirac
63  procedure :: test_randomization => nicas_test_randomization
64  procedure :: test_consistency => nicas_test_consistency
65  procedure :: test_optimality => nicas_test_optimality
66 end type nicas_type
67 
68 private
69 public :: nicas_type
70 
71 contains
72 
73 !----------------------------------------------------------------------
74 ! Subroutine: nicas_alloc
75 ! Purpose: NICAS data allocation
76 !----------------------------------------------------------------------
77 subroutine nicas_alloc(nicas,mpl,nam,bpar,prefix)
78 
79 ! Passed variables
80 class(nicas_type),intent(inout) :: nicas ! NICAS data
81 type(mpl_type),intent(in) :: mpl ! MPI data
82 type(nam_type),intent(in) :: nam ! Namelist
83 type(bpar_type),intent(in) :: bpar ! Block parameters
84 character(len=*),intent(in) :: prefix ! Prefix
85 
86 ! Local variable
87 integer :: ib
88 
89 ! Copy prefix
90 nicas%prefix = prefix
91 
92 ! Allocation
93 allocate(nicas%blk(bpar%nbe))
94 
95 ! Set name
96 do ib=1,bpar%nbe
97  nicas%blk(ib)%ib = ib
98  if (nam%lsqrt) then
99  write(nicas%blk(ib)%name,'(a,i1,a,i4.4,a,i4.4,a,a)') trim(prefix)//'-',nam%mpicom,'-sqrt_',mpl%nproc,'-',mpl%myproc, &
100  & '_',trim(bpar%blockname(ib))
101  else
102  write(nicas%blk(ib)%name,'(a,i1,a,i4.4,a,i4.4,a,a)') trim(prefix)//'-',nam%mpicom,'_',mpl%nproc,'-',mpl%myproc, &
103  & '_',trim(bpar%blockname(ib))
104  end if
105 end do
106 
107 ! Update allocation flag
108 nicas%allocated = .true.
109 
110 end subroutine nicas_alloc
111 
112 !----------------------------------------------------------------------
113 ! Subroutine: nicas_dealloc
114 ! Purpose: NICAS data deallocation
115 !----------------------------------------------------------------------
116 subroutine nicas_dealloc(nicas,nam,geom,bpar)
117 
118 ! Passed variables
119 class(nicas_type),intent(inout) :: nicas ! NICAS data
120 type(nam_type),intent(in) :: nam ! Namelist
121 type(geom_type),intent(in) :: geom ! Geometry
122 type(bpar_type),intent(in) :: bpar ! Block parameters
123 
124 ! Local variable
125 integer :: ib
126 
127 ! Release memory
128 if (allocated(nicas%blk)) then
129  do ib=1,bpar%nbe
130  call nicas%blk(ib)%dealloc(nam,geom)
131  end do
132  deallocate(nicas%blk)
133 end if
134 
135 ! Update allocation flag
136 nicas%allocated = .false.
137 
138 end subroutine nicas_dealloc
139 
140 !----------------------------------------------------------------------
141 ! Subroutine: nicas_read
142 ! Purpose: read NICAS data
143 !----------------------------------------------------------------------
144 subroutine nicas_read(nicas,mpl,nam,geom,bpar)
145 
146 implicit none
147 
148 ! Passed variables
149 class(nicas_type),intent(inout) :: nicas ! NICAS data
150 type(mpl_type),intent(in) :: mpl ! MPI data
151 type(nam_type),intent(in) :: nam ! Namelist
152 type(geom_type),intent(in) :: geom ! Geometry
153 type(bpar_type),intent(in) :: bpar ! Block parameters
154 
155 ! Local variables
156 integer :: ib,il0i,il1,its,il0
157 integer :: info,nl0_test,nc0a_test
158 integer :: ncid,nl0_id,nc0a_id,nc1b_id,nl1_id,nsa_id,nsb_id,nc0d_id,nc0dinv_id
159 integer :: sb_to_c1b_id,sb_to_l1_id
160 integer :: sa_to_sc_id,sb_to_sc_id
161 integer :: norm_id,coef_ens_id
162 character(len=1024) :: filename
163 character(len=1024) :: subr = 'nicas_read'
164 
165 ! Allocation
166 call nicas%alloc(mpl,nam,bpar,'nicas')
167 
168 do ib=1,bpar%nbe
169  if (bpar%B_block(ib)) then
170  ! Open file and get dimensions
171  filename = trim(nam%prefix)//'_'//trim(nicas%blk(ib)%name)//'.nc'
172  call mpl%ncerr(subr,nf90_open(trim(nam%datadir)//'/'//trim(filename),nf90_nowrite,ncid))
173 
174  ! Get dimensions
175  call mpl%ncerr(subr,nf90_inq_dimid(ncid,'nl0',nl0_id))
176  call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nl0_id,len=nl0_test))
177  if (nl0_test/=geom%nl0) call mpl%abort('wrong dimension when reading nicas')
178  if (bpar%nicas_block(ib)) then
179  call mpl%ncerr(subr,nf90_inq_dimid(ncid,'nc0a',nc0a_id))
180  call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nc0a_id,len=nc0a_test))
181  if (nc0a_test/=geom%nc0a) call mpl%abort('wrong dimension when reading nicas')
182  info = nf90_inq_dimid(ncid,'nc1b',nc1b_id)
183  if (info==nf90_noerr) then
184  call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nc1b_id,len=nicas%blk(ib)%nc1b))
185  else
186  nicas%blk(ib)%nc1b = 0
187  end if
188  call mpl%ncerr(subr,nf90_inq_dimid(ncid,'nl1',nl1_id))
189  call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nl1_id,len=nicas%blk(ib)%nl1))
190  info = nf90_inq_dimid(ncid,'nsa',nsa_id)
191  if (info==nf90_noerr) then
192  call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nsa_id,len=nicas%blk(ib)%nsa))
193  else
194  nicas%blk(ib)%nsa = 0
195  end if
196  info = nf90_inq_dimid(ncid,'nsb',nsb_id)
197  if (info==nf90_noerr) then
198  call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nsb_id,len=nicas%blk(ib)%nsb))
199  else
200  nicas%blk(ib)%nsb = 0
201  end if
202  call mpl%ncerr(subr,nf90_get_att(ncid,nf90_global,'nsc',nicas%blk(ib)%nsc))
203  end if
204  if ((ib==bpar%nbe).and.(abs(nam%advmode)==1)) then
205  call mpl%ncerr(subr,nf90_inq_dimid(ncid,'nc0d',nc0d_id))
206  call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nc0d_id,len=nicas%blk(ib)%nc0d))
207  call mpl%ncerr(subr,nf90_inq_dimid(ncid,'nc0dinv',nc0dinv_id))
208  call mpl%ncerr(subr,nf90_inquire_dimension(ncid,nc0dinv_id,len=nicas%blk(ib)%nc0dinv))
209  end if
210 
211  ! Allocation
212  if (bpar%nicas_block(ib)) then
213  if (nicas%blk(ib)%nsb>0) allocate(nicas%blk(ib)%sb_to_c1b(nicas%blk(ib)%nsb))
214  if (nicas%blk(ib)%nsb>0) allocate(nicas%blk(ib)%sb_to_l1(nicas%blk(ib)%nsb))
215  if (nicas%blk(ib)%nsa>0) allocate(nicas%blk(ib)%sa_to_sc(nicas%blk(ib)%nsa))
216  if (nicas%blk(ib)%nsb>0) allocate(nicas%blk(ib)%sb_to_sc(nicas%blk(ib)%nsb))
217  allocate(nicas%blk(ib)%norm(geom%nc0a,geom%nl0))
218  allocate(nicas%blk(ib)%coef_ens(geom%nc0a,geom%nl0))
219  allocate(nicas%blk(ib)%h(geom%nl0i))
220  allocate(nicas%blk(ib)%s(nicas%blk(ib)%nl1))
221  end if
222  if ((ib==bpar%nbe).and.(abs(nam%advmode)==1)) then
223  allocate(nicas%blk(ib)%d(geom%nl0,2:nam%nts))
224  allocate(nicas%blk(ib)%dinv(geom%nl0,2:nam%nts))
225  end if
226 
227  ! Get variable id
228  if (bpar%nicas_block(ib)) then
229  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_inq_varid(ncid,'sb_to_c1b',sb_to_c1b_id))
230  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_inq_varid(ncid,'sb_to_l1',sb_to_l1_id))
231  if (nicas%blk(ib)%nsa>0) call mpl%ncerr(subr,nf90_inq_varid(ncid,'sa_to_sc',sa_to_sc_id))
232  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_inq_varid(ncid,'sb_to_sc',sb_to_sc_id))
233  call mpl%ncerr(subr,nf90_inq_varid(ncid,'norm',norm_id))
234  call mpl%ncerr(subr,nf90_inq_varid(ncid,'coef_ens',coef_ens_id))
235  end if
236 
237  ! Read data
238  if (bpar%nicas_block(ib)) then
239  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_get_var(ncid,sb_to_c1b_id,nicas%blk(ib)%sb_to_c1b))
240  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_get_var(ncid,sb_to_l1_id,nicas%blk(ib)%sb_to_l1))
241  if (nicas%blk(ib)%nsa>0) call mpl%ncerr(subr,nf90_get_var(ncid,sa_to_sc_id,nicas%blk(ib)%sa_to_sc))
242  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_get_var(ncid,sb_to_sc_id,nicas%blk(ib)%sb_to_sc))
243  call mpl%ncerr(subr,nf90_get_var(ncid,norm_id,nicas%blk(ib)%norm))
244  call mpl%ncerr(subr,nf90_get_var(ncid,coef_ens_id,nicas%blk(ib)%coef_ens))
245  call nicas%blk(ib)%com_AB%read(mpl,ncid,'com_AB')
246  call nicas%blk(ib)%com_AC%read(mpl,ncid,'com_AC')
247  nicas%blk(ib)%c%prefix = 'c'
248  call nicas%blk(ib)%c%read(mpl,ncid)
249  do il0i=1,geom%nl0i
250  write(nicas%blk(ib)%h(il0i)%prefix,'(a,i3.3)') 'h_',il0i
251  call nicas%blk(ib)%h(il0i)%read(mpl,ncid)
252  end do
253  nicas%blk(ib)%v%prefix = 'v'
254  call nicas%blk(ib)%v%read(mpl,ncid)
255  do il1=1,nicas%blk(ib)%nl1
256  write(nicas%blk(ib)%s(il1)%prefix,'(a,i3.3)') 's_',il1
257  call nicas%blk(ib)%s(il1)%read(mpl,ncid)
258  end do
259  end if
260  if ((ib==bpar%nbe).and.(abs(nam%advmode)==1)) then
261  call nicas%blk(ib)%com_AD%read(mpl,ncid,'com_AD')
262  call nicas%blk(ib)%com_ADinv%read(mpl,ncid,'com_ADinv')
263  do its=2,nam%nts
264  do il0=1,geom%nl0
265  write(nicas%blk(ib)%d(il0,its)%prefix,'(a,i3.3,a,i2.2)') 'd_',il0,'_',its
266  call nicas%blk(ib)%d(il0,its)%read(mpl,ncid)
267  write(nicas%blk(ib)%dinv(il0,its)%prefix,'(a,i3.3,a,i2.2)') 'dinv_',il0,'_',its
268  call nicas%blk(ib)%dinv(il0,its)%read(mpl,ncid)
269  end do
270  end do
271  end if
272 
273  ! Read main weight
274  call mpl%ncerr(subr,nf90_get_att(ncid,nf90_global,'wgt',nicas%blk(ib)%wgt))
275 
276  ! Close file
277  call mpl%ncerr(subr,nf90_close(ncid))
278  end if
279 end do
280 
281 end subroutine nicas_read
282 
283 !----------------------------------------------------------------------
284 ! Subroutine: nicas_write
285 ! Purpose: write NICAS data
286 !----------------------------------------------------------------------
287 subroutine nicas_write(nicas,mpl,nam,geom,bpar)
288 
289 implicit none
290 
291 ! Passed variables
292 class(nicas_type),intent(in) :: nicas ! NICAS data
293 type(mpl_type),intent(in) :: mpl ! MPI data
294 type(nam_type),intent(in) :: nam ! Namelist
295 type(geom_type),intent(in) :: geom ! Geometry
296 type(bpar_type),intent(in) :: bpar ! Block parameters
297 
298 ! Local variables
299 integer :: ib,il0i,il1,its,il0
300 integer :: ncid,nl0_id,nc0a_id,nc1b_id,nl1_id,nsa_id,nsb_id,nc0d_id,nc0dinv_id
301 integer :: sb_to_c1b_id,sb_to_l1_id
302 integer :: sa_to_sc_id,sb_to_sc_id
303 integer :: norm_id,coef_ens_id
304 character(len=1024) :: filename
305 character(len=1024) :: subr = 'nicas_write'
306 
307 do ib=1,bpar%nbe
308  if (bpar%B_block(ib)) then
309  ! Create file
310  filename = trim(nam%prefix)//'_'//trim(nicas%blk(ib)%name)//'.nc'
311  call mpl%ncerr(subr,nf90_create(trim(nam%datadir)//'/'//trim(filename),or(nf90_clobber,nf90_64bit_offset),ncid))
312 
313  ! Write namelist parameters
314  call nam%ncwrite(mpl,ncid)
315 
316  ! Define dimensions
317  call mpl%ncerr(subr,nf90_def_dim(ncid,'nl0',geom%nl0,nl0_id))
318  if (bpar%nicas_block(ib)) then
319  call mpl%ncerr(subr,nf90_def_dim(ncid,'nc0a',geom%nc0a,nc0a_id))
320  call mpl%ncerr(subr,nf90_put_att(ncid,nf90_global,'nl0i',geom%nl0i))
321  if (nicas%blk(ib)%nc1b>0) call mpl%ncerr(subr,nf90_def_dim(ncid,'nc1b',nicas%blk(ib)%nc1b,nc1b_id))
322  call mpl%ncerr(subr,nf90_def_dim(ncid,'nl1',nicas%blk(ib)%nl1,nl1_id))
323  if (nicas%blk(ib)%nsa>0) call mpl%ncerr(subr,nf90_def_dim(ncid,'nsa',nicas%blk(ib)%nsa,nsa_id))
324  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_def_dim(ncid,'nsb',nicas%blk(ib)%nsb,nsb_id))
325  call mpl%ncerr(subr,nf90_put_att(ncid,nf90_global,'nsc',nicas%blk(ib)%nsc))
326  end if
327  if ((ib==bpar%nbe).and.nam%displ_diag) then
328  call mpl%ncerr(subr,nf90_def_dim(ncid,'nc0d',nicas%blk(ib)%nc0d,nc0d_id))
329  call mpl%ncerr(subr,nf90_def_dim(ncid,'nc0dinv',nicas%blk(ib)%nc0dinv,nc0dinv_id))
330  end if
331 
332  ! Write main weight
333  call mpl%ncerr(subr,nf90_put_att(ncid,nf90_global,'wgt',nicas%blk(ib)%wgt))
334 
335  ! Define variables
336  if (bpar%nicas_block(ib)) then
337  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_def_var(ncid,'sb_to_c1b',nf90_int,(/nsb_id/),sb_to_c1b_id))
338  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_def_var(ncid,'sb_to_l1',nf90_int,(/nsb_id/),sb_to_l1_id))
339  if (nicas%blk(ib)%nsa>0) call mpl%ncerr(subr,nf90_def_var(ncid,'sa_to_sc',nf90_int,(/nsa_id/),sa_to_sc_id))
340  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_def_var(ncid,'sb_to_sc',nf90_int,(/nsb_id/),sb_to_sc_id))
341  call mpl%ncerr(subr,nf90_def_var(ncid,'norm',ncfloat,(/nc0a_id,nl0_id/),norm_id))
342  call mpl%ncerr(subr,nf90_def_var(ncid,'coef_ens',ncfloat,(/nc0a_id,nl0_id/),coef_ens_id))
343 
344  if (nicas%blk(ib)%nsa>0) call mpl%ncerr(subr,nf90_put_att(ncid,sa_to_sc_id,'_FillValue',msvali))
345  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_put_att(ncid,sb_to_sc_id,'_FillValue',msvali))
346  call mpl%ncerr(subr,nf90_put_att(ncid,norm_id,'_FillValue',msvalr))
347  call mpl%ncerr(subr,nf90_put_att(ncid,coef_ens_id,'_FillValue',msvalr))
348  end if
349 
350  ! End definition mode
351  call mpl%ncerr(subr,nf90_enddef(ncid))
352 
353  ! Write variables
354  if (bpar%nicas_block(ib)) then
355  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_put_var(ncid,sb_to_c1b_id,nicas%blk(ib)%sb_to_c1b))
356  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_put_var(ncid,sb_to_l1_id,nicas%blk(ib)%sb_to_l1))
357  if (nicas%blk(ib)%nsa>0) call mpl%ncerr(subr,nf90_put_var(ncid,sa_to_sc_id,nicas%blk(ib)%sa_to_sc))
358  if (nicas%blk(ib)%nsb>0) call mpl%ncerr(subr,nf90_put_var(ncid,sb_to_sc_id,nicas%blk(ib)%sb_to_sc))
359  call mpl%ncerr(subr,nf90_put_var(ncid,norm_id,nicas%blk(ib)%norm))
360  call mpl%ncerr(subr,nf90_put_var(ncid,coef_ens_id,nicas%blk(ib)%coef_ens))
361  call nicas%blk(ib)%com_AB%write(mpl,ncid)
362  call nicas%blk(ib)%com_AC%write(mpl,ncid)
363  call nicas%blk(ib)%c%write(mpl,ncid)
364  do il0i=1,geom%nl0i
365  call nicas%blk(ib)%h(il0i)%write(mpl,ncid)
366  end do
367  call nicas%blk(ib)%v%write(mpl,ncid)
368  do il1=1,nicas%blk(ib)%nl1
369  call nicas%blk(ib)%s(il1)%write(mpl,ncid)
370  end do
371  end if
372  if ((ib==bpar%nbe).and.nam%displ_diag) then
373  call nicas%blk(ib)%com_AD%write(mpl,ncid)
374  call nicas%blk(ib)%com_ADinv%write(mpl,ncid)
375  do its=2,nam%nts
376  do il0=1,geom%nl0
377  call nicas%blk(ib)%d(il0,its)%write(mpl,ncid)
378  call nicas%blk(ib)%dinv(il0,its)%write(mpl,ncid)
379  end do
380  end do
381  end if
382 
383  ! Close file
384  call mpl%ncerr(subr,nf90_close(ncid))
385  end if
386 end do
387 
388 end subroutine nicas_write
389 
390 !----------------------------------------------------------------------
391 ! Subroutine: nicas_write_mpi_summary
392 ! Purpose: write NICAS MPI related data summary
393 !----------------------------------------------------------------------
394 subroutine nicas_write_mpi_summary(nicas,mpl,nam,geom,bpar)
395 
396 implicit none
397 
398 ! Passed variables
399 class(nicas_type),intent(in) :: nicas ! NICAS data
400 type(mpl_type),intent(in) :: mpl ! MPI data
401 type(nam_type),intent(in) :: nam ! Namelist
402 type(geom_type),intent(in) :: geom ! Geometry
403 type(bpar_type),intent(in) :: bpar ! Block parameters
404 
405 ! Local variables
406 integer :: ib,is,ic1,il1
407 integer :: ncid,nc0_id,nc1_id,nl1_id
408 integer :: lon_id,lat_id,c0_to_proc_id,c1_to_c0_id,l1_to_l0_id,lcheck_id
409 real(kind_real),allocatable :: lcheck(:,:)
410 character(len=1024) :: filename
411 character(len=1024) :: subr = 'nicas_write_mpi_summary'
412 
413 do ib=1,bpar%nbe
414  if (bpar%nicas_block(ib)) then
415  ! Allocation
416  allocate(lcheck(nicas%blk(ib)%nc1,nicas%blk(ib)%nl1))
417 
418  ! Create summary file
419  filename = trim(nam%prefix)//'_'//trim(nicas%blk(ib)%name)//'_summary.nc'
420  call mpl%ncerr(subr,nf90_create(trim(nam%datadir)//'/'//trim(filename),or(nf90_clobber,nf90_64bit_offset),ncid))
421 
422  ! Write namelist parameters
423  call nam%ncwrite(mpl,ncid)
424 
425  ! Define dimensions
426  call mpl%ncerr(subr,nf90_def_dim(ncid,'nc0',geom%nc0,nc0_id))
427  call mpl%ncerr(subr,nf90_def_dim(ncid,'nc1',nicas%blk(ib)%nc1,nc1_id))
428  call mpl%ncerr(subr,nf90_def_dim(ncid,'nl1',nicas%blk(ib)%nl1,nl1_id))
429 
430  ! Define variables
431  call mpl%ncerr(subr,nf90_def_var(ncid,'lon',ncfloat,(/nc0_id/),lon_id))
432  call mpl%ncerr(subr,nf90_def_var(ncid,'lat',ncfloat,(/nc0_id/),lat_id))
433  call mpl%ncerr(subr,nf90_def_var(ncid,'c0_to_proc',nf90_int,(/nc0_id/),c0_to_proc_id))
434  call mpl%ncerr(subr,nf90_def_var(ncid,'c1_to_c0',nf90_int,(/nc1_id/),c1_to_c0_id))
435  call mpl%ncerr(subr,nf90_def_var(ncid,'l1_to_l0',nf90_int,(/nl1_id/),l1_to_l0_id))
436  call mpl%ncerr(subr,nf90_def_var(ncid,'lcheck',ncfloat,(/nc1_id,nl1_id/),lcheck_id))
437 
438  ! End definition mode
439  call mpl%ncerr(subr,nf90_enddef(ncid))
440 
441  ! Write variables
442  call mpl%ncerr(subr,nf90_put_var(ncid,lon_id,geom%lon*rad2deg))
443  call mpl%ncerr(subr,nf90_put_var(ncid,lat_id,geom%lat*rad2deg))
444  call mpl%ncerr(subr,nf90_put_var(ncid,c0_to_proc_id,geom%c0_to_proc))
445  call mpl%ncerr(subr,nf90_put_var(ncid,c1_to_c0_id,nicas%blk(ib)%c1_to_c0))
446  call mpl%ncerr(subr,nf90_put_var(ncid,l1_to_l0_id,nicas%blk(ib)%l1_to_l0))
447  call msr(lcheck)
448  do is=1,nicas%blk(ib)%ns
449  ic1 = nicas%blk(ib)%s_to_c1(is)
450  il1 = nicas%blk(ib)%s_to_l1(is)
451  if (nicas%blk(ib)%lcheck_sa(is)) then
452  lcheck(ic1,il1) = 1.0
453  elseif (nicas%blk(ib)%lcheck_sb(is)) then
454  lcheck(ic1,il1) = 2.0
455  elseif (nicas%blk(ib)%lcheck_sc(is)) then
456  lcheck(ic1,il1) = 3.0
457  else
458  lcheck(ic1,il1) = 4.0
459  end if
460  end do
461  call mpl%ncerr(subr,nf90_put_var(ncid,lcheck_id,lcheck))
462 
463  ! Close summary file
464  call mpl%ncerr(subr,nf90_close(ncid))
465 
466  ! Release memory
467  deallocate(lcheck)
468  end if
469 end do
470 
471 end subroutine nicas_write_mpi_summary
472 
473 !----------------------------------------------------------------------
474 ! Subroutine: nicas_run_nicas
475 ! Purpose: NICAS driver
476 !----------------------------------------------------------------------
477 subroutine nicas_run_nicas(nicas,mpl,rng,nam,geom,bpar,cmat)
478 
479 implicit none
480 
481 ! Passed variables
482 class(nicas_type),intent(inout) :: nicas ! NICAS data
483 type(mpl_type),intent(inout) :: mpl ! MPI data
484 type(rng_type),intent(inout) :: rng ! Random number generator
485 type(nam_type),intent(inout) :: nam ! Namelist
486 type(geom_type),intent(inout) :: geom ! Geometry
487 type(bpar_type),intent(in) :: bpar ! Block parameters
488 type(cmat_type),intent(in) :: cmat ! C matrix data
489 
490 ! Local variables
491 integer :: ib
492 
493 ! Allocation
494 call nicas%alloc(mpl,nam,bpar,'nicas')
495 
496 ! Compute NICAS parameters
497 write(mpl%info,'(a)') '-------------------------------------------------------------------'
498 write(mpl%info,'(a)') '--- Compute NICAS parameters'
499 call flush(mpl%info)
500 
501 do ib=1,bpar%nbe
502  if (bpar%nicas_block(ib).or.((ib==bpar%nbe).and.nam%displ_diag)) then
503  write(mpl%info,'(a)') '-------------------------------------------------------------------'
504  write(mpl%info,'(a)') '--- Block: '//trim(bpar%blockname(ib))
505  call flush(mpl%info)
506  end if
507 
508  ! NICAS parameters
509  if (bpar%nicas_block(ib)) call nicas%blk(ib)%compute_parameters(mpl,rng,nam,geom,cmat%blk(ib))
510 
511  ! Advection
512  if ((ib==bpar%nbe).and.nam%displ_diag) call nicas%blk(ib)%compute_adv(mpl,rng,nam,geom,cmat%blk(ib))
513 
514  ! Coefficient
515  if (bpar%B_block(ib)) then
516  ! Copy weights
517  nicas%blk(ib)%wgt = cmat%blk(ib)%wgt
518  if (bpar%nicas_block(ib)) then
519  allocate(nicas%blk(ib)%coef_ens(geom%nc0a,geom%nl0))
520  nicas%blk(ib)%coef_ens = cmat%blk(ib)%coef_ens
521  end if
522  end if
523 end do
524 
525 ! Write NICAS parameters
526 write(mpl%info,'(a)') '-------------------------------------------------------------------'
527 write(mpl%info,'(a)') '--- Write NICAS parameters'
528 call flush(mpl%info)
529 call nicas%write(mpl,nam,geom,bpar)
530 
531 ! Write NICAS MPI summary
532 write(mpl%info,'(a)') '-------------------------------------------------------------------'
533 write(mpl%info,'(a)') '--- Write NICAS MPI summary'
534 call flush(mpl%info)
535 call nicas%write_mpi_summary(mpl,nam,geom,bpar)
536 
537 end subroutine nicas_run_nicas
538 
539 !----------------------------------------------------------------------
540 ! Subroutine: nicas_run_nicas_tests
541 ! Purpose: NICAS tests driver
542 !----------------------------------------------------------------------
543 subroutine nicas_run_nicas_tests(nicas,mpl,rng,nam,geom,bpar,io,cmat,ens)
544 
545 implicit none
546 
547 ! Passed variables
548 class(nicas_type),intent(inout) :: nicas ! NICAS data
549 type(mpl_type),intent(inout) :: mpl ! MPI data
550 type(rng_type),intent(inout) :: rng ! Random number generator
551 type(nam_type),intent(inout) :: nam ! Namelist
552 type(geom_type),intent(inout) :: geom ! Geometry
553 type(bpar_type),intent(in) :: bpar ! Block parameters
554 type(io_type),intent(in) :: io ! I/O
555 type(cmat_type),intent(in) :: cmat ! C matrix data
556 type(ens_type),intent(in) :: ens ! Ensemble
557 
558 ! Local variables
559 integer :: ib
560 
561 if (nam%check_adjoints) then
562  ! Test adjoint
563  write(mpl%info,'(a)') '-------------------------------------------------------------------'
564  write(mpl%info,'(a)') '--- Test NICAS adjoint'
565  call flush(mpl%info)
566 
567  do ib=1,bpar%nbe
568  if (bpar%nicas_block(ib)) then
569  write(mpl%info,'(a)') '-------------------------------------------------------------------'
570  write(mpl%info,'(a)') '--- Block: '//trim(bpar%blockname(ib))
571  call flush(mpl%info)
572  call nicas%blk(ib)%test_adjoint(mpl,rng,nam,geom)
573  end if
574  end do
575 
576  ! Test NICAS adjoint
577  write(mpl%info,'(a)') '-------------------------------------------------------------------'
578  write(mpl%info,'(a)') '--- Test NICAS adjoint'
579  call flush(mpl%info)
580  call nicas%test_adjoint(mpl,rng,nam,geom,bpar,ens)
581 end if
582 
583 if (nam%check_pos_def) then
584  ! Test NICAS positive definiteness
585  write(mpl%info,'(a)') '-------------------------------------------------------------------'
586  write(mpl%info,'(a)') '--- Test NICAS positive definiteness'
587  call flush(mpl%info)
588 
589  do ib=1,bpar%nbe
590  if (bpar%nicas_block(ib)) then
591  write(mpl%info,'(a)') '-------------------------------------------------------------------'
592  write(mpl%info,'(a)') '--- Block: '//trim(bpar%blockname(ib))
593  call flush(mpl%info)
594  call nicas%blk(ib)%test_pos_def(mpl,rng,nam,geom)
595  end if
596  end do
597 end if
598 
599 if (nam%check_sqrt) then
600  ! Test NICAS full/square-root equivalence
601  write(mpl%info,'(a)') '-------------------------------------------------------------------'
602  write(mpl%info,'(a)') '--- Test NICAS full/square-root equivalence'
603  call flush(mpl%info)
604 
605  do ib=1,bpar%nbe
606  if (bpar%nicas_block(ib)) then
607  write(mpl%info,'(a)') '-------------------------------------------------------------------'
608  write(mpl%info,'(a)') '--- Block: '//trim(bpar%blockname(ib))
609  call flush(mpl%info)
610  call nicas%blk(ib)%test_sqrt(mpl,rng,nam,geom,bpar,io,cmat%blk(ib))
611  end if
612  end do
613 
614  ! Test NICAS full/square-root equivalence
615  write(mpl%info,'(a)') '-------------------------------------------------------------------'
616  write(mpl%info,'(a)') '--- Test NICAS full/square-root equivalence'
617  call flush(mpl%info)
618  call nicas%test_sqrt(mpl,rng,nam,geom,bpar,io,cmat,ens)
619 end if
620 
621 if (nam%check_dirac) then
622  ! Apply NICAS to diracs
623  write(mpl%info,'(a)') '-------------------------------------------------------------------'
624  write(mpl%info,'(a)') '--- Apply NICAS to diracs'
625  call flush(mpl%info)
626 
627  do ib=1,bpar%nbe
628  if (bpar%nicas_block(ib)) then
629  write(mpl%info,'(a)') '-------------------------------------------------------------------'
630  write(mpl%info,'(a)') '--- Block: '//trim(bpar%blockname(ib))
631  call flush(mpl%info)
632  call nicas%blk(ib)%test_dirac(mpl,nam,geom,bpar,io)
633  end if
634  end do
635 
636  ! Apply NICAS to diracs
637  write(mpl%info,'(a)') '-------------------------------------------------------------------'
638  write(mpl%info,'(a)') '--- Apply NICAS to diracs'
639  call flush(mpl%info)
640  call nicas%test_dirac(mpl,nam,geom,bpar,io,ens)
641 end if
642 
643 if (nam%check_randomization) then
644  ! Test NICAS randomization
645  write(mpl%info,'(a)') '-------------------------------------------------------------------'
646  write(mpl%info,'(a)') '--- Test NICAS randomization'
647  call flush(mpl%info)
648  call nicas%test_randomization(mpl,rng,nam,geom,bpar,io)
649 end if
650 
651 if (nam%check_consistency) then
652  ! Test HDIAG-NICAS consistency
653  write(mpl%info,'(a)') '-------------------------------------------------------------------'
654  write(mpl%info,'(a)') '--- Test HDIAG-NICAS consistency'
655  call flush(mpl%info)
656  call nicas%test_consistency(mpl,rng,nam,geom,bpar,io,cmat)
657 end if
658 
659 if (nam%check_optimality) then
660  ! Test HDIAG optimality
661  write(mpl%info,'(a)') '-------------------------------------------------------------------'
662  write(mpl%info,'(a)') '--- Test HDIAG optimality'
663  call flush(mpl%info)
664  call nicas%test_optimality(mpl,rng,nam,geom,bpar,io)
665 end if
666 
667 end subroutine nicas_run_nicas_tests
668 
669 !----------------------------------------------------------------------
670 ! Subroutine: nicas_alloc_cv
671 ! Purpose: control vector allocation
672 !----------------------------------------------------------------------
673 subroutine nicas_alloc_cv(nicas,bpar,cv,getsizeonly)
674 
675 implicit none
676 
677 ! Passed variables
678 class(nicas_type),intent(in) :: nicas ! NICAS data
679 type(bpar_type),intent(in) :: bpar ! Block parameters
680 type(cv_type),intent(inout) :: cv ! Control vector
681 logical,intent(in),optional :: getsizeonly ! Flag to get the control variable size only (no allocation)
682 
683 ! Local variables
684 integer :: ib
685 logical :: lgetsizeonly
686 
687 ! Check flag existence
688 lgetsizeonly = .false.
689 if (present(getsizeonly)) lgetsizeonly = getsizeonly
690 
691 ! Allocation
692 allocate(cv%blk(bpar%nbe))
693 
694 ! Initialization
695 cv%n = 0
696 cv%nbe = bpar%nbe
697 
698 do ib=1,bpar%nbe
699  if (bpar%cv_block(ib)) then
700  ! Copy block size
701  cv%blk(ib)%n = nicas%blk(ib)%nsa
702 
703  ! Update total size
704  cv%n = cv%n+nicas%blk(ib)%nsa
705 
706  if (.not.lgetsizeonly) then
707  ! Allocation
708  allocate(cv%blk(ib)%alpha(cv%blk(ib)%n))
709 
710  ! Initialization
711  call msr(cv%blk(ib)%alpha)
712  end if
713  else
714  ! Set zero size
715  cv%blk(ib)%n = 0
716  end if
717 end do
718 
719 end subroutine nicas_alloc_cv
720 
721 !----------------------------------------------------------------------
722 ! Subroutine: nicas_random_cv
723 ! Purpose: generate a random control vector
724 !----------------------------------------------------------------------
725 subroutine nicas_random_cv(nicas,rng,bpar,cv)
726 
727 implicit none
728 
729 ! Passed variables
730 class(nicas_type),intent(in) :: nicas ! NICAS data
731 type(rng_type),intent(inout) :: rng ! Random number generator
732 type(bpar_type),intent(in) :: bpar ! Block parameters
733 type(cv_type),intent(out) :: cv ! Control vector
734 
735 ! Local variables
736 integer :: ib
737 
738 ! Allocation
739 call nicas%alloc_cv(bpar,cv)
740 
741 ! Random initialization
742 do ib=1,bpar%nbe
743  if (bpar%cv_block(ib)) call rng%rand_gau(cv%blk(ib)%alpha)
744 end do
745 
746 end subroutine nicas_random_cv
747 
748 !----------------------------------------------------------------------
749 ! Subroutine: nicas_apply
750 ! Purpose: apply NICAS
751 !----------------------------------------------------------------------
752 subroutine nicas_apply(nicas,mpl,nam,geom,bpar,fld)
753 
754 implicit none
755 
756 ! Passed variables
757 class(nicas_type),intent(in) :: nicas ! NICAS data
758 type(mpl_type),intent(in) :: mpl ! MPI data
759 type(nam_type),target,intent(in) :: nam ! Namelist
760 type(geom_type),target,intent(in) :: geom ! Geometry
761 type(bpar_type),target,intent(in) :: bpar ! Block parameters
762 real(kind_real),intent(inout) :: fld(geom%nc0a,geom%nl0,nam%nv,nam%nts) ! Field
763 
764 ! Local variable
765 integer :: ib,its,iv,jv,il0,ic0a
766 real(kind_real) :: prod,prod_tot
767 real(kind_real),allocatable :: fld_3d(:,:),fld_4d(:,:,:),fld_4d_tmp(:,:,:)
768 real(kind_real),allocatable :: wgt(:,:),wgt_diag(:)
769 real(kind_real),allocatable :: fld_save(:,:,:,:)
770 
771 if (pos_def_test) then
772  ! Save field for positive-definiteness test
773  allocate(fld_save(geom%nc0a,geom%nl0,nam%nv,nam%nts))
774  fld_save = fld
775 end if
776 
777 ! Adjoint advection
778 if (nam%advmode==1) call nicas%blk(bpar%nbe)%apply_adv_ad(mpl,nam,geom,fld)
779 
780 select case (nam%strategy)
781 case ('common')
782  ! Allocation
783  allocate(fld_3d(geom%nc0a,geom%nl0))
784 
785  ! Sum product over variables and timeslots
786  fld_3d = 0.0
787  !$omp parallel do schedule(static) private(il0,ic0a,its,iv)
788  do il0=1,geom%nl0
789  do ic0a=1,geom%nc0a
790  do its=1,nam%nts
791  do iv=1,nam%nv
792  fld_3d(ic0a,il0) = fld_3d(ic0a,il0)+fld(ic0a,il0,iv,its)
793  end do
794  end do
795  end do
796  end do
797  !$omp end parallel do
798 
799  ! Apply common ensemble coefficient square-root
800  !$omp parallel do schedule(static) private(il0,ic0a)
801  do il0=1,geom%nl0
802  do ic0a=1,geom%nc0a
803  if (geom%mask_c0a(ic0a,il0)) fld_3d(ic0a,il0) = fld_3d(ic0a,il0)*sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
804  end do
805  end do
806  !$omp end parallel do
807 
808  ! Apply common NICAS
809  call nicas%blk(bpar%nbe)%apply(mpl,nam,geom,fld_3d)
810 
811  ! Apply common ensemble coefficient square-root
812  !$omp parallel do schedule(static) private(il0,ic0a)
813  do il0=1,geom%nl0
814  do ic0a=1,geom%nc0a
815  if (geom%mask_c0a(ic0a,il0)) fld_3d(ic0a,il0) = fld_3d(ic0a,il0)*sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
816  end do
817  end do
818  !$omp end parallel do
819 
820  ! Build final vector
821  !$omp parallel do schedule(static) private(il0,ic0a,its,iv)
822  do il0=1,geom%nl0
823  do ic0a=1,geom%nc0a
824  do its=1,nam%nts
825  do iv=1,nam%nv
826  fld(ic0a,il0,iv,its) = fld_3d(ic0a,il0)
827  end do
828  end do
829  end do
830  end do
831  !$omp end parallel do
832 case ('common_univariate')
833  ! Allocation
834  allocate(fld_4d(geom%nc0a,geom%nl0,nam%nv))
835 
836  ! Sum product over timeslots
837  fld_4d = 0.0
838  do its=1,nam%nts
839  fld_4d = fld_4d+fld(:,:,:,its)
840  end do
841 
842  do iv=1,nam%nv
843  ! Apply common ensemble coefficient square-root
844  !$omp parallel do schedule(static) private(il0,ic0a)
845  do il0=1,geom%nl0
846  do ic0a=1,geom%nc0a
847  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
848  & *sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
849  end do
850  end do
851  !$omp end parallel do
852 
853  ! Apply common NICAS
854  call nicas%blk(bpar%nbe)%apply(mpl,nam,geom,fld_4d(:,:,iv))
855 
856  ! Apply common ensemble coefficient square-root
857  !$omp parallel do schedule(static) private(il0,ic0a)
858  do il0=1,geom%nl0
859  do ic0a=1,geom%nc0a
860  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
861  & *sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
862  end do
863  end do
864  !$omp end parallel do
865  end do
866 
867  ! Build final vector
868  do its=1,nam%nts
869  fld(:,:,:,its) = fld_4d
870  end do
871 case ('common_weighted')
872  ! Allocation
873  allocate(fld_4d(geom%nc0a,geom%nl0,nam%nv))
874  allocate(fld_4d_tmp(geom%nc0a,geom%nl0,nam%nv))
875  allocate(wgt(nam%nv,nam%nv))
876  allocate(wgt_diag(nam%nv))
877 
878  ! Copy weights
879  do ib=1,bpar%nb
880  if (bpar%B_block(ib)) then
881  ! Variable indices
882  iv = bpar%b_to_v1(ib)
883  jv = bpar%b_to_v2(ib)
884  if (isnotmsr(nicas%blk(ib)%wgt)) then
885  wgt(iv,jv) = nicas%blk(ib)%wgt
886  else
887  wgt(iv,jv) = 0.0
888  end if
889  if (iv==jv) wgt_diag(iv) = wgt(iv,iv)
890  end if
891  end do
892 
893  ! Normalize weights
894  do iv=1,nam%nv
895  do jv=1,nam%nv
896  wgt(iv,jv) = wgt(iv,jv)/sqrt(wgt_diag(iv)*wgt_diag(jv))
897  end do
898  end do
899 
900  ! Sum product over timeslots
901  fld_4d = 0.0
902  do its=1,nam%nts
903  fld_4d = fld_4d+fld(:,:,:,its)
904  end do
905 
906  do iv=1,nam%nv
907  ! Apply common ensemble coefficient square-root
908  !$omp parallel do schedule(static) private(il0,ic0a)
909  do il0=1,geom%nl0
910  do ic0a=1,geom%nc0a
911  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
912  & *sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
913  end do
914  end do
915  !$omp end parallel do
916 
917  ! Apply common NICAS
918  call nicas%blk(bpar%nbe)%apply(mpl,nam,geom,fld_4d(:,:,iv))
919 
920  ! Apply common ensemble coefficient square-root
921  !$omp parallel do schedule(static) private(il0,ic0a)
922  do il0=1,geom%nl0
923  do ic0a=1,geom%nc0a
924  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
925  & *sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
926  end do
927  end do
928  !$omp end parallel do
929  end do
930 
931  ! Apply weights
932  fld_4d_tmp = 0.0
933  do iv=1,nam%nv
934  do jv=1,nam%nv
935  fld_4d_tmp(:,:,iv) = fld_4d_tmp(:,:,iv)+wgt(iv,jv)*fld_4d(:,:,jv)
936  end do
937  end do
938 
939  ! Build final vector
940  do its=1,nam%nts
941  fld(:,:,:,its) = fld_4d_tmp
942  end do
943 case ('specific_univariate')
944  ! Allocation
945  allocate(fld_4d(geom%nc0a,geom%nl0,nam%nv))
946 
947  ! Sum product over timeslots
948  fld_4d = 0.0
949  do its=1,nam%nts
950  fld_4d = fld_4d+fld(:,:,:,its)
951  end do
952 
953  do ib=1,bpar%nb
954  if (bpar%nicas_block(ib)) then
955  ! Variable index
956  iv = bpar%b_to_v1(ib)
957 
958  ! Apply common ensemble coefficient square-root
959  !$omp parallel do schedule(static) private(il0,ic0a)
960  do il0=1,geom%nl0
961  do ic0a=1,geom%nc0a
962  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
963  & *sqrt(nicas%blk(ib)%coef_ens(ic0a,il0))
964  end do
965  end do
966  !$omp end parallel do
967 
968  ! Apply specific NICAS (same for all timeslots)
969  call nicas%blk(ib)%apply(mpl,nam,geom,fld_4d(:,:,iv))
970 
971  ! Apply common ensemble coefficient square-root
972  !$omp parallel do schedule(static) private(il0,ic0a)
973  do il0=1,geom%nl0
974  do ic0a=1,geom%nc0a
975  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
976  & *sqrt(nicas%blk(ib)%coef_ens(ic0a,il0))
977  end do
978  end do
979  !$omp end parallel do
980  end if
981  end do
982 
983  ! Build final vector
984  do its=1,nam%nts
985  fld(:,:,:,its) = fld_4d
986  end do
987 case ('specific_multivariate')
988  call mpl%abort('specific multivariate strategy should not be called from apply_NICAS (lsqrt required)')
989 end select
990 
991 if (pos_def_test) then
992  ! Positive-definiteness test
993  prod = sum(fld_save*fld)
994  call mpl%f_comm%allreduce(prod,prod_tot,fckit_mpi_sum())
995  if (prod_tot<0.0) call mpl%abort('negative result in nicas_apply')
996 end if
997 
998 ! Advection
999 if (nam%advmode==1) call nicas%blk(bpar%nbe)%apply_adv(mpl,nam,geom,fld)
1000 
1001 end subroutine nicas_apply
1002 
1003 !----------------------------------------------------------------------
1004 ! Subroutine: nicas_apply_from_sqrt
1005 ! Purpose: apply NICAS from square-root
1006 !----------------------------------------------------------------------
1007 subroutine nicas_apply_from_sqrt(nicas,mpl,nam,geom,bpar,fld)
1008 
1009 implicit none
1010 
1011 ! Passed variables
1012 class(nicas_type),intent(in) :: nicas ! NICAS data
1013 type(mpl_type),intent(in) :: mpl ! MPI data
1014 type(nam_type),target,intent(in) :: nam ! Namelist
1015 type(geom_type),target,intent(in) :: geom ! Geometry
1016 type(bpar_type),target,intent(in) :: bpar ! Block parameters
1017 real(kind_real),intent(inout) :: fld(geom%nc0a,geom%nl0,nam%nv,nam%nts) ! Field
1018 
1019 ! Local variable
1020 real(kind_real) :: prod,prod_tot
1021 real(kind_real),allocatable :: fld_save(:,:,:,:)
1022 type(cv_type) :: cv
1023 
1024 if (pos_def_test) then
1025  ! Save field for positivity test
1026  allocate(fld_save(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1027  fld_save = fld
1028 end if
1029 
1030 ! Apply square-root adjoint
1031 call nicas%apply_sqrt_ad(mpl,nam,geom,bpar,fld,cv)
1032 
1033 ! Apply square-root
1034 call nicas%apply_sqrt(mpl,nam,geom,bpar,cv,fld)
1035 
1036 if (pos_def_test) then
1037  ! Positivity test
1038  prod = sum(fld_save*fld)
1039  call mpl%f_comm%allreduce(prod,prod_tot,fckit_mpi_sum())
1040  if (prod_tot<0.0) call mpl%abort('negative result in nicas_apply')
1041 end if
1042 
1043 end subroutine nicas_apply_from_sqrt
1044 
1045 !----------------------------------------------------------------------
1046 ! Subroutine: nicas_apply_sqrt
1047 ! Purpose: apply NICAS square-root
1048 !----------------------------------------------------------------------
1049 subroutine nicas_apply_sqrt(nicas,mpl,nam,geom,bpar,cv,fld)
1050 
1051 implicit none
1052 
1053 ! Passed variables
1054 class(nicas_type),intent(in) :: nicas ! NICAS data
1055 type(mpl_type),intent(in) :: mpl ! MPI data
1056 type(nam_type),target,intent(in) :: nam ! Namelist
1057 type(geom_type),target,intent(in) :: geom ! Geometry
1058 type(bpar_type),target,intent(in) :: bpar ! Block parameters
1059 type(cv_type),intent(in) :: cv ! Control variable
1060 real(kind_real),intent(out) :: fld(geom%nc0a,geom%nl0,nam%nv,nam%nts) ! Field
1061 
1062 ! Local variable
1063 integer :: ib,its,iv,jv,ic0a,il0
1064 real(kind_real),allocatable :: fld_3d(:,:),fld_4d(:,:,:),fld_4d_tmp(:,:,:)
1065 real(kind_real),allocatable :: wgt(:,:),wgt_diag(:),wgt_u(:,:)
1066 
1067 select case (nam%strategy)
1068 case ('common')
1069  ! Allocation
1070  allocate(fld_3d(geom%nc0a,geom%nl0))
1071 
1072  ! Apply common NICAS
1073  call nicas%blk(bpar%nbe)%apply_sqrt(mpl,geom,cv%blk(bpar%nbe)%alpha,fld_3d)
1074 
1075  ! Apply common ensemble coefficient square-root
1076  !$omp parallel do schedule(static) private(il0,ic0a)
1077  do il0=1,geom%nl0
1078  do ic0a=1,geom%nc0a
1079  if (geom%mask_c0a(ic0a,il0)) fld_3d(ic0a,il0) = fld_3d(ic0a,il0)*sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
1080  end do
1081  end do
1082  !$omp end parallel do
1083 
1084  ! Build final vector
1085  do its=1,nam%nts
1086  do iv=1,nam%nv
1087  fld(:,:,iv,its) = fld_3d
1088  end do
1089  end do
1090 case ('common_univariate')
1091  ! Allocation
1092  allocate(fld_4d(geom%nc0a,geom%nl0,nam%nv))
1093 
1094  do ib=1,bpar%nb
1095  if (bpar%cv_block(ib)) then
1096  ! Variable index
1097  iv = bpar%b_to_v1(ib)
1098 
1099  ! Apply specific NICAS (same for all timeslots)
1100  call nicas%blk(bpar%nbe)%apply_sqrt(mpl,geom,cv%blk(ib)%alpha,fld_4d(:,:,iv))
1101 
1102  ! Apply common ensemble coefficient square-root
1103  !$omp parallel do schedule(static) private(il0,ic0a)
1104  do il0=1,geom%nl0
1105  do ic0a=1,geom%nc0a
1106  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
1107  & *sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
1108  end do
1109  end do
1110  !$omp end parallel do
1111  end if
1112  end do
1113 
1114  ! Build final vector
1115  do its=1,nam%nts
1116  fld(:,:,:,its) = fld_4d
1117  end do
1118 case ('common_weighted')
1119  ! Allocation
1120  allocate(fld_4d(geom%nc0a,geom%nl0,nam%nv))
1121  allocate(fld_4d_tmp(geom%nc0a,geom%nl0,nam%nv))
1122  allocate(wgt(nam%nv,nam%nv))
1123  allocate(wgt_diag(nam%nv))
1124  allocate(wgt_u(nam%nv,nam%nv))
1125 
1126  ! Copy weights
1127  do ib=1,bpar%nb
1128  if (bpar%B_block(ib)) then
1129  ! Variable indices
1130  iv = bpar%b_to_v1(ib)
1131  jv = bpar%b_to_v2(ib)
1132  if (isnotmsr(nicas%blk(ib)%wgt)) then
1133  wgt(iv,jv) = nicas%blk(ib)%wgt
1134  else
1135  wgt(iv,jv) = 0.0
1136  end if
1137  if (iv==jv) wgt_diag(iv) = wgt(iv,iv)
1138  end if
1139  end do
1140 
1141  ! Normalize weights
1142  do iv=1,nam%nv
1143  do jv=1,nam%nv
1144  wgt(iv,jv) = wgt(iv,jv)/sqrt(wgt_diag(iv)*wgt_diag(jv))
1145  end do
1146  end do
1147 
1148  ! Cholesky decomposition
1149  call cholesky(mpl,nam%nv,wgt,wgt_u)
1150 
1151  do ib=1,bpar%nb
1152  if (bpar%cv_block(ib)) then
1153  ! Variable index
1154  iv = bpar%b_to_v1(ib)
1155 
1156  ! Apply specific NICAS (same for all timeslots)
1157  call nicas%blk(bpar%nbe)%apply_sqrt(mpl,geom,cv%blk(ib)%alpha,fld_4d(:,:,iv))
1158 
1159  ! Apply common ensemble coefficient square-root
1160  !$omp parallel do schedule(static) private(il0,ic0a)
1161  do il0=1,geom%nl0
1162  do ic0a=1,geom%nc0a
1163  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
1164  & *sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
1165  end do
1166  end do
1167  !$omp end parallel do
1168  end if
1169  end do
1170 
1171  ! Apply weights
1172  fld_4d_tmp = 0.0
1173  do iv=1,nam%nv
1174  do jv=1,iv
1175  fld_4d_tmp(:,:,iv) = fld_4d_tmp(:,:,iv)+wgt_u(iv,jv)*fld_4d(:,:,jv)
1176  end do
1177  end do
1178 
1179  ! Build final vector
1180  do its=1,nam%nts
1181  fld(:,:,:,its) = fld_4d_tmp
1182  end do
1183 case ('specific_univariate')
1184  ! Allocation
1185  allocate(fld_4d(geom%nc0a,geom%nl0,nam%nv))
1186 
1187  do ib=1,bpar%nb
1188  if (bpar%nicas_block(ib)) then
1189  ! Variable index
1190  iv = bpar%b_to_v1(ib)
1191 
1192  ! Apply specific NICAS (same for all timeslots)
1193  call nicas%blk(ib)%apply_sqrt(mpl,geom,cv%blk(ib)%alpha,fld_4d(:,:,iv))
1194 
1195  ! Apply specific ensemble coefficient square-root
1196  !$omp parallel do schedule(static) private(il0,ic0a)
1197  do il0=1,geom%nl0
1198  do ic0a=1,geom%nc0a
1199  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
1200  & *sqrt(nicas%blk(ib)%coef_ens(ic0a,il0))
1201  end do
1202  end do
1203  !$omp end parallel do
1204  end if
1205  end do
1206 
1207  ! Build final vector
1208  do its=1,nam%nts
1209  fld(:,:,:,its) = fld_4d
1210  end do
1211 case ('specific_multivariate')
1212  ! Allocation
1213  allocate(fld_4d(geom%nc0a,geom%nl0,nam%nv))
1214 
1215  do ib=1,bpar%nb
1216  if (bpar%nicas_block(ib)) then
1217  ! Variable index
1218  iv = bpar%b_to_v1(ib)
1219 
1220  ! Apply specific NICAS (same for all timeslots)
1221  call nicas%blk(ib)%apply_sqrt(mpl,geom,cv%blk(bpar%nbe)%alpha,fld_4d(:,:,iv))
1222 
1223  ! Apply specific ensemble coefficient square-root
1224  !$omp parallel do schedule(static) private(il0,ic0a)
1225  do il0=1,geom%nl0
1226  do ic0a=1,geom%nc0a
1227  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
1228  & *sqrt(nicas%blk(ib)%coef_ens(ic0a,il0))
1229  end do
1230  end do
1231  !$omp end parallel do
1232  end if
1233  end do
1234 
1235  ! Build final vector
1236  do its=1,nam%nts
1237  fld(:,:,:,its) = fld_4d
1238  end do
1239 end select
1240 
1241 ! Advection
1242 if (nam%advmode==1) call nicas%blk(bpar%nbe)%apply_adv(mpl,nam,geom,fld)
1243 
1244 end subroutine nicas_apply_sqrt
1245 
1246 !----------------------------------------------------------------------
1247 ! Subroutine: nicas_apply_sqrt_ad
1248 ! Purpose: apply NICAS square-root, adjoint
1249 !----------------------------------------------------------------------
1250 subroutine nicas_apply_sqrt_ad(nicas,mpl,nam,geom,bpar,fld,cv)
1251 
1252 implicit none
1253 
1254 ! Passed variables
1255 class(nicas_type),intent(in) :: nicas ! NICAS data
1256 type(mpl_type),intent(in) :: mpl ! MPI data
1257 type(nam_type),target,intent(in) :: nam ! Namelist
1258 type(geom_type),target,intent(in) :: geom ! Geometry
1259 type(bpar_type),target,intent(in) :: bpar ! Block parameters
1260 real(kind_real),intent(in) :: fld(geom%nc0a,geom%nl0,nam%nv,nam%nts) ! Field
1261 type(cv_type),intent(out) :: cv ! Control variable
1262 
1263 ! Local variable
1264 integer :: ib,its,iv,jv,ic0a,il0
1265 real(kind_real),allocatable :: fld_3d(:,:),fld_4d(:,:,:),fld_4d_tmp(:,:,:),fld_5d(:,:,:,:)
1266 real(kind_real),allocatable :: wgt(:,:),wgt_diag(:),wgt_u(:,:)
1267 type(cv_type) :: cv_tmp
1268 
1269 ! Allocation
1270 allocate(fld_5d(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1271 call nicas%alloc_cv(bpar,cv)
1272 
1273 ! Copy
1274 fld_5d = fld
1275 
1276 ! Adjoint advection
1277 if (nam%advmode==1) call nicas%blk(bpar%nbe)%apply_adv_ad(mpl,nam,geom,fld_5d)
1278 
1279 select case (nam%strategy)
1280 case ('common')
1281  ! Allocation
1282  allocate(fld_3d(geom%nc0a,geom%nl0))
1283 
1284  ! Sum product over variables and timeslots
1285  fld_3d = 0.0
1286  do its=1,nam%nts
1287  do iv=1,nam%nv
1288  fld_3d = fld_3d+fld_5d(:,:,iv,its)
1289  end do
1290  end do
1291 
1292  ! Apply common ensemble coefficient square-root
1293  !$omp parallel do schedule(static) private(il0,ic0a)
1294  do il0=1,geom%nl0
1295  do ic0a=1,geom%nc0a
1296  if (geom%mask_c0a(ic0a,il0)) fld_3d(ic0a,il0) = fld_3d(ic0a,il0)*sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
1297  end do
1298  end do
1299  !$omp end parallel do
1300 
1301  ! Apply common NICAS
1302  call nicas%blk(bpar%nbe)%apply_sqrt_ad(mpl,geom,fld_3d,cv%blk(bpar%nbe)%alpha)
1303 case ('common_univariate')
1304  ! Allocation
1305  allocate(fld_4d(geom%nc0a,geom%nl0,nam%nv))
1306 
1307  ! Sum product over timeslots
1308  fld_4d = 0.0
1309  do its=1,nam%nts
1310  fld_4d = fld_4d+fld_5d(:,:,:,its)
1311  end do
1312 
1313  do ib=1,bpar%nb
1314  if (bpar%cv_block(ib)) then
1315  ! Variable index
1316  iv = bpar%b_to_v1(ib)
1317 
1318  ! Apply common ensemble coefficient square-root
1319  !$omp parallel do schedule(static) private(il0,ic0a)
1320  do il0=1,geom%nl0
1321  do ic0a=1,geom%nc0a
1322  if (geom%mask_c0a(ic0a,il0)) fld_4d_tmp(ic0a,il0,iv) = fld_4d_tmp(ic0a,il0,iv) &
1323  & *sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
1324  end do
1325  end do
1326  !$omp end parallel do
1327 
1328  ! Apply specific NICAS (same for all timeslots)
1329  call nicas%blk(bpar%nbe)%apply_sqrt_ad(mpl,geom,fld_4d_tmp(:,:,iv),cv%blk(ib)%alpha)
1330  end if
1331  end do
1332 case ('common_weighted')
1333  ! Allocation
1334  allocate(fld_4d(geom%nc0a,geom%nl0,nam%nv))
1335  allocate(fld_4d_tmp(geom%nc0a,geom%nl0,nam%nv))
1336  allocate(wgt(nam%nv,nam%nv))
1337  allocate(wgt_diag(nam%nv))
1338  allocate(wgt_u(nam%nv,nam%nv))
1339 
1340  ! Copy weights
1341  do ib=1,bpar%nb
1342  if (bpar%B_block(ib)) then
1343  ! Variable indices
1344  iv = bpar%b_to_v1(ib)
1345  jv = bpar%b_to_v2(ib)
1346  if (isnotmsr(nicas%blk(ib)%wgt)) then
1347  wgt(iv,jv) = nicas%blk(ib)%wgt
1348  else
1349  wgt(iv,jv) = 0.0
1350  end if
1351  if (iv==jv) wgt_diag(iv) = wgt(iv,iv)
1352  end if
1353  end do
1354 
1355  ! Normalize weights
1356  do iv=1,nam%nv
1357  do jv=1,nam%nv
1358  wgt(iv,jv) = wgt(iv,jv)/sqrt(wgt_diag(iv)*wgt_diag(jv))
1359  end do
1360  end do
1361 
1362  ! Cholesky decomposition
1363  call cholesky(mpl,nam%nv,wgt,wgt_u)
1364 
1365  ! Sum product over timeslots
1366  fld_4d = 0.0
1367  do its=1,nam%nts
1368  fld_4d = fld_4d+fld_5d(:,:,:,its)
1369  end do
1370 
1371  ! Apply weights
1372  fld_4d_tmp = 0.0
1373  do iv=1,nam%nv
1374  do jv=iv,nam%nv
1375  fld_4d_tmp(:,:,iv) = fld_4d_tmp(:,:,iv)+wgt_u(iv,jv)*fld_4d(:,:,jv)
1376  end do
1377  end do
1378 
1379  do ib=1,bpar%nb
1380  if (bpar%cv_block(ib)) then
1381  ! Variable index
1382  iv = bpar%b_to_v1(ib)
1383 
1384  ! Apply common ensemble coefficient square-root
1385  !$omp parallel do schedule(static) private(il0,ic0a)
1386  do il0=1,geom%nl0
1387  do ic0a=1,geom%nc0a
1388  if (geom%mask_c0a(ic0a,il0)) fld_4d_tmp(ic0a,il0,iv) = fld_4d_tmp(ic0a,il0,iv) &
1389  & *sqrt(nicas%blk(bpar%nbe)%coef_ens(ic0a,il0))
1390  end do
1391  end do
1392  !$omp end parallel do
1393 
1394  ! Apply specific NICAS (same for all timeslots)
1395  call nicas%blk(bpar%nbe)%apply_sqrt_ad(mpl,geom,fld_4d_tmp(:,:,iv),cv%blk(ib)%alpha)
1396  end if
1397  end do
1398 case ('specific_univariate')
1399  ! Allocation
1400  allocate(fld_4d(geom%nc0a,geom%nl0,nam%nv))
1401 
1402  ! Sum product over timeslots
1403  fld_4d = 0.0
1404  do its=1,nam%nts
1405  fld_4d = fld_4d+fld_5d(:,:,:,its)
1406  end do
1407 
1408  do ib=1,bpar%nb
1409  if (bpar%nicas_block(ib)) then
1410  ! Variable index
1411  iv = bpar%b_to_v1(ib)
1412 
1413  ! Apply common ensemble coefficient square-root
1414  !$omp parallel do schedule(static) private(il0,ic0a)
1415  do il0=1,geom%nl0
1416  do ic0a=1,geom%nc0a
1417  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
1418  & *sqrt(nicas%blk(ib)%coef_ens(ic0a,il0))
1419  end do
1420  end do
1421  !$omp end parallel do
1422 
1423  ! Apply specific NICAS (same for all timeslots)
1424  call nicas%blk(ib)%apply_sqrt_ad(mpl,geom,fld_4d(:,:,iv),cv%blk(ib)%alpha)
1425  end if
1426  end do
1427 case ('specific_multivariate')
1428  ! Allocation
1429  allocate(fld_4d(geom%nc0a,geom%nl0,nam%nv))
1430  call nicas%alloc_cv(bpar,cv_tmp)
1431 
1432  ! Initialization
1433  cv%blk(bpar%nbe)%alpha = 0.0
1434 
1435  ! Sum product over timeslots
1436  fld_4d = 0.0
1437  do its=1,nam%nts
1438  fld_4d = fld_4d+fld_5d(:,:,:,its)
1439  end do
1440 
1441  do ib=1,bpar%nb
1442  if (bpar%nicas_block(ib)) then
1443  ! Variable index
1444  iv = bpar%b_to_v1(ib)
1445 
1446  ! Apply common ensemble coefficient square-root
1447  !$omp parallel do schedule(static) private(il0,ic0a)
1448  do il0=1,geom%nl0
1449  do ic0a=1,geom%nc0a
1450  if (geom%mask_c0a(ic0a,il0)) fld_4d(ic0a,il0,iv) = fld_4d(ic0a,il0,iv) &
1451  & *sqrt(nicas%blk(ib)%coef_ens(ic0a,il0))
1452  end do
1453  end do
1454  !$omp end parallel do
1455 
1456  ! Apply specific NICAS (same for all timeslots)
1457  call nicas%blk(ib)%apply_sqrt_ad(mpl,geom,fld_4d(:,:,iv),cv_tmp%blk(bpar%nbe)%alpha)
1458 
1459  ! Sum control variable
1460  cv%blk(bpar%nbe)%alpha = cv%blk(bpar%nbe)%alpha+cv_tmp%blk(bpar%nbe)%alpha
1461  end if
1462  end do
1463 end select
1464 
1465 end subroutine nicas_apply_sqrt_ad
1466 
1467 !----------------------------------------------------------------------
1468 ! Subroutine: nicas_randomize
1469 ! Purpose: randomize NICAS from square-root
1470 !----------------------------------------------------------------------
1471 subroutine nicas_randomize(nicas,mpl,rng,nam,geom,bpar,ne,ens)
1472 
1473 implicit none
1474 
1475 ! Passed variables
1476 class(nicas_type),intent(in) :: nicas ! NICAS data
1477 type(mpl_type),intent(in) :: mpl ! MPI data
1478 type(rng_type),intent(inout) :: rng ! Random number generator
1479 type(nam_type),target,intent(in) :: nam ! Namelist
1480 type(geom_type),target,intent(in) :: geom ! Geometry
1481 type(bpar_type),target,intent(in) :: bpar ! Blocal parameters
1482 integer,intent(in) :: ne ! Number of members
1483 type(ens_type),intent(out) :: ens ! Ensemble
1484 
1485 ! Local variable
1486 integer :: ie,ic0a,il0,its,iv
1487 real(kind_real) :: norm,mean(geom%nc0a,geom%nl0,nam%nv,nam%nts),std(geom%nc0a,geom%nl0,nam%nv,nam%nts)
1488 type(cv_type) :: cv_ens(ne)
1489 
1490 ! Allocation
1491 call ens%alloc(nam,geom,ne,1)
1492 
1493 do ie=1,ne
1494  ! Generate random control vector
1495  call nicas%random_cv(rng,bpar,cv_ens(ie))
1496 
1497  ! Apply square-root
1498  call nicas%apply_sqrt(mpl,nam,geom,bpar,cv_ens(ie),ens%fld(:,:,:,:,ie))
1499 end do
1500 
1501 ! Remove mean
1502 mean = sum(ens%fld,dim=5)/real(ne,kind_real)
1503 do ie=1,ne
1504  ens%fld(:,:,:,:,ie) = ens%fld(:,:,:,:,ie)-mean
1505 end do
1506 
1507 ! Compute standard deviation
1508 norm = real(ne-1,kind_real)
1509 !$omp parallel do schedule(static) private(its,iv,il0,ic0a)
1510 do its=1,nam%nts
1511  do iv=1,nam%nv
1512  do il0=1,geom%nl0
1513  do ic0a=1,geom%nc0a
1514  if (geom%mask_c0a(ic0a,il0)) std(ic0a,il0,iv,its) = sqrt(sum(ens%fld(ic0a,il0,iv,its,:)**2)) &
1515  & /norm
1516  end do
1517  end do
1518  end do
1519 end do
1520 !$omp end parallel do
1521 
1522 ! Normalize perturbations
1523 !$omp parallel do schedule(static) private(its,iv,il0,ic0a)
1524 do its=1,nam%nts
1525  do iv=1,nam%nv
1526  do il0=1,geom%nl0
1527  do ic0a=1,geom%nc0a
1528  if (geom%mask_c0a(ic0a,il0)) ens%fld(ic0a,il0,iv,its,:) = ens%fld(ic0a,il0,iv,its,:) &
1529  & /std(ic0a,il0,iv,its)
1530  end do
1531  end do
1532  end do
1533 end do
1534 !$omp end parallel do
1535 
1536 end subroutine nicas_randomize
1537 
1538 !----------------------------------------------------------------------
1539 ! Subroutine: nicas_apply_bens
1540 ! Purpose: apply localized ensemble covariance
1541 !----------------------------------------------------------------------
1542 subroutine nicas_apply_bens(nicas,mpl,nam,geom,bpar,ens,fld)
1543 
1544 implicit none
1545 
1546 ! Passed variables
1547 class(nicas_type),intent(in) :: nicas ! NICAS data
1548 type(mpl_type),intent(in) :: mpl ! MPI data
1549 type(nam_type),target,intent(in) :: nam ! Namelist
1550 type(geom_type),target,intent(in) :: geom ! Geometry
1551 type(bpar_type),target,intent(in) :: bpar ! Blocal parameters
1552 type(ens_type),intent(in) :: ens ! Ensemble
1553 real(kind_real),intent(inout) :: fld(geom%nc0a,geom%nl0,nam%nv,nam%nts) ! Field
1554 
1555 ! Local variable
1556 integer :: ie
1557 real(kind_real) :: fld_copy(geom%nc0a,geom%nl0,nam%nv,nam%nts),fld_tmp(geom%nc0a,geom%nl0,nam%nv,nam%nts)
1558 real(kind_real) :: pert(geom%nc0a,geom%nl0,nam%nv,nam%nts)
1559 
1560 ! Copy field
1561 fld_copy = fld
1562 
1563 ! Adjoint advection
1564 if (nam%advmode==-1) call nicas%blk(bpar%nbe)%apply_adv_ad(mpl,nam,geom,fld_copy)
1565 
1566 ! Apply localized ensemble covariance formula
1567 fld = 0.0
1568 do ie=1,nam%ens1_ne
1569  ! Compute perturbation
1570  pert = ens%fld(:,:,:,:,ie)
1571 
1572  ! Inverse advection
1573  if (nam%advmode==-1) call nicas%blk(bpar%nbe)%apply_adv_inv(mpl,nam,geom,pert)
1574 
1575  ! Schur product
1576  fld_tmp = pert*fld_copy
1577 
1578  ! Apply NICAS
1579  if (nam%lsqrt) then
1580  call nicas%apply_from_sqrt(mpl,nam,geom,bpar,fld_tmp)
1581  else
1582  call nicas%apply(mpl,nam,geom,bpar,fld_tmp)
1583  end if
1584 
1585  ! Schur product
1586  fld = fld+fld_tmp*pert
1587 
1588  ! Normalization
1589  fld = fld/real(nam%ens1_ne-1,kind_real)
1590 end do
1591 
1592 ! Advection
1593 if (nam%advmode==-1) call nicas%blk(bpar%nbe)%apply_adv(mpl,nam,geom,fld)
1594 
1595 end subroutine nicas_apply_bens
1596 
1597 !----------------------------------------------------------------------
1598 ! Subroutine: nicas_apply_bens_noloc
1599 ! Purpose: apply ensemble covariance, without localization
1600 !----------------------------------------------------------------------
1601 subroutine nicas_apply_bens_noloc(nicas,mpl,nam,geom,ens,fld)
1602 
1603 implicit none
1604 
1605 ! Passed variables
1606 class(nicas_type),intent(in) :: nicas ! NICAS data
1607 type(mpl_type),intent(in) :: mpl ! MPI data
1608 type(nam_type),target,intent(in) :: nam ! Namelist
1609 type(geom_type),target,intent(in) :: geom ! Geometry
1610 type(ens_type),intent(in) :: ens ! Ensemble
1611 real(kind_real),intent(inout) :: fld(geom%nc0a,geom%nl0,nam%nv,nam%nts) ! Field
1612 
1613 ! Local variable
1614 integer :: ie,ic0a,il0,iv,its
1615 real(kind_real) :: alpha,norm
1616 real(kind_real) :: fld_copy(geom%nc0a,geom%nl0,nam%nv,nam%nts)
1617 real(kind_real) :: pert(geom%nc0a,geom%nl0,nam%nv,nam%nts)
1618 character(len=1024) :: dum
1619 
1620 ! Initialization
1621 fld_copy = fld
1622 
1623 ! Apply localized ensemble covariance formula
1624 fld = 0.0
1625 norm = sqrt(real(nam%ens1_ne-1,kind_real))
1626 do ie=1,nam%ens1_ne
1627  ! Compute perturbation
1628  !$omp parallel do schedule(static) private(its,iv,il0,ic0a)
1629  do its=1,nam%nts
1630  do iv=1,nam%nv
1631  do il0=1,geom%nl0
1632  do ic0a=1,geom%nc0a
1633  if (geom%mask_c0a(ic0a,il0)) pert(ic0a,il0,iv,its) = ens%fld(ic0a,il0,iv,its,ie)/norm
1634  end do
1635  end do
1636  end do
1637  end do
1638  !$omp end parallel do
1639 
1640  ! Dot product
1641  call mpl%dot_prod(pert,fld_copy,alpha)
1642 
1643  ! Schur product
1644  fld = fld+alpha*pert
1645 end do
1646 
1647 ! To avoid compilation warnings
1648 dum = nicas%prefix
1649 
1650 end subroutine nicas_apply_bens_noloc
1651 
1652 !----------------------------------------------------------------------
1653 ! Subroutine: nicas_test_adjoint
1654 ! Purpose: test NICAS adjoint
1655 !----------------------------------------------------------------------
1656 subroutine nicas_test_adjoint(nicas,mpl,rng,nam,geom,bpar,ens)
1657 
1658 implicit none
1659 
1660 ! Passed variables
1661 class(nicas_type),intent(in) :: nicas ! NICAS data
1662 type(mpl_type),intent(in) :: mpl ! MPI data
1663 type(rng_type),intent(inout) :: rng ! Random number generator
1664 type(nam_type),intent(in) :: nam ! Namelist
1665 type(geom_type),intent(in) :: geom ! Geometry
1666 type(bpar_type),intent(in) :: bpar ! Block parameters
1667 type(ens_type),intent(in) :: ens ! Ensemble
1668 
1669 ! Local variables
1670 real(kind_real) :: sum1,sum2
1671 real(kind_real),allocatable :: fld1_loc(:,:,:,:),fld1_adv(:,:,:,:),fld1_bens(:,:,:,:),fld1_save(:,:,:,:)
1672 real(kind_real),allocatable :: fld2_loc(:,:,:,:),fld2_adv(:,:,:,:),fld2_bens(:,:,:,:),fld2_save(:,:,:,:)
1673 
1674 ! Allocation
1675 allocate(fld1_save(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1676 allocate(fld2_save(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1677 allocate(fld1_loc(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1678 allocate(fld2_loc(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1679 if (ens%ne>0) then
1680  allocate(fld1_bens(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1681  allocate(fld2_bens(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1682 end if
1683 
1684 ! Generate random field
1685 call rng%rand_real(0.0_kind_real,1.0_kind_real,fld1_save)
1686 call rng%rand_real(0.0_kind_real,1.0_kind_real,fld2_save)
1687 
1688 ! Adjoint test
1689 fld1_loc = fld1_save
1690 fld2_loc = fld2_save
1691 if (nam%lsqrt) then
1692  call nicas%apply_from_sqrt(mpl,nam,geom,bpar,fld1_loc)
1693  call nicas%apply_from_sqrt(mpl,nam,geom,bpar,fld2_loc)
1694 else
1695  call nicas%apply(mpl,nam,geom,bpar,fld1_loc)
1696  call nicas%apply(mpl,nam,geom,bpar,fld2_loc)
1697 end if
1698 if (abs(nam%advmode)==1) then
1699  fld1_adv = fld1_save
1700  fld2_adv = fld2_save
1701  call nicas%blk(bpar%nbe)%apply_adv(mpl,nam,geom,fld1_adv)
1702  call nicas%blk(bpar%nbe)%apply_adv_ad(mpl,nam,geom,fld2_adv)
1703 end if
1704 if (ens%ne>0) then
1705  fld1_bens = fld1_save
1706  fld2_bens = fld2_save
1707  call nicas%apply_bens(mpl,nam,geom,bpar,ens,fld1_bens)
1708  call nicas%apply_bens(mpl,nam,geom,bpar,ens,fld2_bens)
1709 end if
1710 
1711 ! Print result
1712 call mpl%dot_prod(fld1_loc,fld2_save,sum1)
1713 call mpl%dot_prod(fld2_loc,fld1_save,sum2)
1714 write(mpl%info,'(a7,a,e15.8,a,e15.8,a,e15.8)') '','NICAS adjoint test: ', &
1715  & sum1,' / ',sum2,' / ',2.0*abs(sum1-sum2)/abs(sum1+sum2)
1716 call flush(mpl%info)
1717 if (abs(nam%advmode)==1) then
1718  call mpl%dot_prod(fld1_adv,fld2_save,sum1)
1719  call mpl%dot_prod(fld2_adv,fld1_save,sum2)
1720  write(mpl%info,'(a7,a,e15.8,a,e15.8,a,e15.8)') '','Advection adjoint test: ', &
1721  & sum1,' / ',sum2,' / ',2.0*abs(sum1-sum2)/abs(sum1+sum2)
1722  call flush(mpl%info)
1723 end if
1724 if (ens%ne>0) then
1725  call mpl%dot_prod(fld1_bens,fld2_save,sum1)
1726  call mpl%dot_prod(fld2_bens,fld1_save,sum2)
1727  write(mpl%info,'(a7,a,e15.8,a,e15.8,a,e15.8)') '','Ensemble B adjoint test: ', &
1728  & sum1,' / ',sum2,' / ',2.0*abs(sum1-sum2)/abs(sum1+sum2)
1729  call flush(mpl%info)
1730 end if
1731 
1732 end subroutine nicas_test_adjoint
1733 
1734 !----------------------------------------------------------------------
1735 ! Subroutine: nicas_test_sqrt
1736 ! Purpose: test full/square-root equivalence
1737 !----------------------------------------------------------------------
1738 subroutine nicas_test_sqrt(nicas,mpl,rng,nam,geom,bpar,io,cmat,ens)
1739 
1740 implicit none
1741 
1742 ! Passed variables
1743 class(nicas_type),intent(in) :: nicas ! NICAS data
1744 type(mpl_type),intent(inout) :: mpl ! MPI data
1745 type(rng_type),intent(inout) :: rng ! Random number generator
1746 type(nam_type),intent(inout),target :: nam ! Namelist
1747 type(geom_type),intent(in),target :: geom ! Geometry
1748 type(bpar_type),intent(in) :: bpar ! Block parameters
1749 type(io_type),intent(in) :: io ! I/O
1750 type(cmat_type),intent(in) :: cmat ! C matrix data
1751 type(ens_type),intent(in) :: ens ! Ensemble
1752 
1753 ! Local variables
1754 integer :: ib,iv
1755 real(kind_real),allocatable :: fld_loc(:,:,:,:),fld_loc_sqrt(:,:,:,:)
1756 real(kind_real),allocatable :: fld_bens(:,:,:,:),fld_bens_sqrt(:,:,:,:)
1757 character(len=1024) :: varname(nam%nv)
1758 type(nicas_type) :: nicas_other
1759 
1760 ! Allocation
1761 allocate(fld_loc(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1762 allocate(fld_loc_sqrt(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1763 if (ens%ne>0) then
1764  allocate(fld_bens(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1765  allocate(fld_bens_sqrt(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1766 end if
1767 
1768 ! Generate random field
1769 call rng%rand_real(-1.0_kind_real,1.0_kind_real,fld_loc)
1770 fld_loc_sqrt = fld_loc
1771 if (ens%ne>0) then
1772  call rng%rand_real(-1.0_kind_real,1.0_kind_real,fld_bens)
1773  fld_bens_sqrt = fld_bens
1774 end if
1775 
1776 ! Apply NICAS, initial version
1777 if (nam%lsqrt) then
1778  call nicas%apply_from_sqrt(mpl,nam,geom,bpar,fld_loc_sqrt)
1779 else
1780  call nicas%apply(mpl,nam,geom,bpar,fld_loc)
1781 end if
1782 if (ens%ne>0) then
1783  if (nam%lsqrt) then
1784  call nicas%apply_bens(mpl,nam,geom,bpar,ens,fld_bens_sqrt)
1785  else
1786  call nicas%apply_bens(mpl,nam,geom,bpar,ens,fld_bens)
1787  end if
1788 end if
1789 
1790 ! Switch lsqrt
1791 nam%lsqrt = .not.nam%lsqrt
1792 
1793 ! Allocation
1794 call nicas_other%alloc(mpl,nam,bpar,'nicas_other')
1795 
1796 ! Prepare nicas, other version
1797 do ib=1,bpar%nbe
1798  if (bpar%nicas_block(ib)) then
1799  ! Compute NICAS parameters
1800  call nicas_other%blk(ib)%compute_parameters(mpl,rng,nam,geom,cmat%blk(ib))
1801  end if
1802 
1803  if (bpar%B_block(ib)) then
1804  ! Copy weights
1805  nicas_other%blk(ib)%wgt = cmat%blk(ib)%wgt
1806  if (bpar%nicas_block(ib)) then
1807  allocate(nicas_other%blk(ib)%coef_ens(geom%nc0a,geom%nl0))
1808  nicas_other%blk(ib)%coef_ens = cmat%blk(ib)%coef_ens
1809  end if
1810  end if
1811 end do
1812 
1813 ! Apply NICAS, other version
1814 if (nam%lsqrt) then
1815  call nicas_other%apply_from_sqrt(mpl,nam,geom,bpar,fld_loc_sqrt)
1816 else
1817  call nicas_other%apply(mpl,nam,geom,bpar,fld_loc)
1818 end if
1819 if (ens%ne>0) then
1820  if (nam%lsqrt) then
1821  call nicas_other%apply_bens(mpl,nam,geom,bpar,ens,fld_bens_sqrt)
1822  else
1823  call nicas_other%apply_bens(mpl,nam,geom,bpar,ens,fld_bens)
1824  end if
1825 end if
1826 
1827 ! Compute dirac
1828 do iv=1,nam%nv
1829  varname(iv) = nam%varname(iv)
1830  nam%varname(iv) = trim(varname(iv))//'_sqrt'
1831 end do
1832 if (nam%check_dirac) call nicas_other%test_dirac(mpl,nam,geom,bpar,io,ens)
1833 do iv=1,nam%nv
1834  nam%varname(iv) = varname(iv)
1835 end do
1836 
1837 ! Reset lsqrt value
1838 nam%lsqrt = .not.nam%lsqrt
1839 
1840 ! Print difference
1841 write(mpl%info,'(a7,a,f6.1,a)') '','NICAS full / square-root error : ', &
1842  & sqrt(sum((fld_loc_sqrt-fld_loc)**2)/sum(fld_loc**2))*100.0,'%'
1843 if (ens%ne>0) write(mpl%info,'(a7,a,f6.1,a)') '','Ensemble B full / square-root error: ', &
1844  & sqrt(sum((fld_bens_sqrt-fld_bens)**2)/sum(fld_bens**2))*100.0,'%'
1845 call flush(mpl%info)
1846 
1847 end subroutine nicas_test_sqrt
1848 
1849 !----------------------------------------------------------------------
1850 ! Subroutine: nicas_test_dirac
1851 ! Purpose: apply NICAS to diracs
1852 !----------------------------------------------------------------------
1853 subroutine nicas_test_dirac(nicas,mpl,nam,geom,bpar,io,ens)
1854 
1855 implicit none
1856 
1857 ! Passed variables
1858 class(nicas_type),intent(in) :: nicas ! NICAS data
1859 type(mpl_type),intent(inout) :: mpl ! MPI data
1860 type(nam_type),intent(in) :: nam ! Namelist
1861 type(geom_type),intent(in) :: geom ! Geometry
1862 type(bpar_type),intent(in) :: bpar ! Block parameters
1863 type(io_type),intent(in) :: io ! I/O
1864 type(ens_type),intent(in) :: ens ! Ensemble
1865 
1866 ! Local variables
1867 integer :: idir,iv,its
1868 real(kind_real),allocatable :: fld(:,:,:,:),fld_loc(:,:,:,:),fld_bens(:,:,:,:)
1869 character(len=2) :: itschar
1870 character(len=1024) :: filename
1871 
1872 ! Allocation
1873 allocate(fld(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1874 
1875 ! Generate dirac field
1876 fld = 0.0
1877 do idir=1,geom%ndir
1878  if (geom%iprocdir(idir)==mpl%myproc) fld(geom%ic0adir(idir),geom%il0dir(idir),geom%ivdir(idir),geom%itsdir(idir)) = 1.0
1879 end do
1880 
1881 ! Allocation
1882 allocate(fld_loc(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1883 if (ens%ne>0) allocate(fld_bens(geom%nc0a,geom%nl0,nam%nv,nam%nts))
1884 
1885 ! Apply NICAS to dirac
1886 fld_loc = fld
1887 if (nam%lsqrt) then
1888  call nicas%apply_from_sqrt(mpl,nam,geom,bpar,fld_loc)
1889 else
1890  call nicas%apply(mpl,nam,geom,bpar,fld_loc)
1891 end if
1892 
1893 if ((ens%ne>0).and.(trim(nam%method)/='cor')) then
1894  ! Apply localized ensemble covariance
1895  fld_bens = fld
1896  call nicas%apply_bens(mpl,nam,geom,bpar,ens,fld_bens)
1897 end if
1898 
1899 ! Write field
1900 filename = trim(nam%prefix)//'_dirac'
1901 do its=1,nam%nts
1902  write(itschar,'(i2.2)') its
1903  do iv=1,nam%nv
1904  call io%fld_write(mpl,nam,geom,filename,trim(nam%varname(iv))//'_'//itschar,fld_loc(:,:,iv,its))
1905  if ((ens%ne>0).and.(trim(nam%method)/='cor')) call io%fld_write(mpl,nam,geom,filename, &
1906  & trim(nam%varname(iv))//'_'//itschar//'_Bens',fld_bens(:,:,iv,its))
1907  end do
1908 end do
1909 
1910 end subroutine nicas_test_dirac
1911 
1912 !----------------------------------------------------------------------
1913 ! Subroutine: nicas_test_randomization
1914 ! Purpose: test NICAS randomization method with respect to theoretical error statistics
1915 !----------------------------------------------------------------------
1916 subroutine nicas_test_randomization(nicas,mpl,rng,nam,geom,bpar,io)
1917 
1918 implicit none
1919 
1920 ! Passed variables
1921 class(nicas_type),intent(in) :: nicas ! NICAS data
1922 type(mpl_type),intent(inout) :: mpl ! MPI data
1923 type(rng_type),intent(inout) :: rng ! Random number generator
1924 type(nam_type),intent(inout) :: nam ! Namelist variables
1925 type(geom_type),intent(in) :: geom ! Geometry
1926 type(bpar_type),intent(in) :: bpar ! Block parameters
1927 type(io_type),intent(in) :: io ! I/O
1928 
1929 ! Local variables
1930 integer :: ifac,itest,nefac(nfac),ens1_ne,iv,its
1931 real(kind_real) :: fld_ref(geom%nc0a,geom%nl0,nam%nv,nam%nts,ntest),fld_save(geom%nc0a,geom%nl0,nam%nv,nam%nts,ntest)
1932 real(kind_real) :: fld(geom%nc0a,geom%nl0,nam%nv,nam%nts),mse(ntest,nfac),mse_th(ntest,nfac)
1933 character(len=2) :: itschar
1934 character(len=4) :: nechar,itestchar
1935 character(len=1024) :: filename
1936 type(ens_type) :: ens
1937 
1938 ! Define test vectors
1939 write(mpl%info,'(a4,a)') '','Define test vectors'
1940 call flush(mpl%info)
1941 call define_test_vectors(mpl,rng,nam,geom,ntest,fld_save)
1942 
1943 ! Apply NICAS to test vectors
1944 write(mpl%info,'(a4,a)') '','Apply NICAS to test vectors'
1945 call flush(mpl%info)
1946 fld_ref = fld_save
1947 do itest=1,ntest
1948  call nicas%apply_from_sqrt(mpl,nam,geom,bpar,fld_ref(:,:,:,:,itest))
1949 end do
1950 
1951 ! Write first 10 test vectors
1952 do itest=1,min(ntest,10)
1953  ! Write field
1954  write(itestchar,'(i4.4)') itest
1955  filename = trim(nam%prefix)//'_randomize_'//itestchar
1956  do its=1,nam%nts
1957  write(itschar,'(i2.2)') its
1958  do iv=1,nam%nv
1959  call io%fld_write(mpl,nam,geom,filename,trim(nam%varname(iv))//'_ref_'//itschar,fld_ref(:,:,iv,its,itest))
1960  end do
1961  end do
1962 end do
1963 
1964 ! Save namelist variables
1965 ens1_ne = nam%ens1_ne
1966 
1967 write(mpl%info,'(a4,a)') '','Test randomization for various ensemble sizes:'
1968 call flush(mpl%info)
1969 do ifac=1,nfac
1970  ! Ensemble size
1971  nefac(ifac) = max(int(5.0*real(ifac,kind_real)/real(nfac,kind_real)*real(ne_rand,kind_real)),3)
1972  nam%ens1_ne = nefac(ifac)
1973  write(nechar,'(i4.4)') nefac(ifac)
1974 
1975  ! Randomize ensemble
1976  call nicas%randomize(mpl,rng,nam,geom,bpar,nefac(ifac),ens)
1977 
1978  do itest=1,ntest
1979  ! Test NICAS
1980  fld = fld_save(:,:,:,:,itest)
1981  call nicas%apply_bens_noloc(mpl,nam,geom,ens,fld)
1982 
1983  ! RMSE
1984  mse(itest,ifac) = sum((fld-fld_ref(:,:,:,:,itest))**2)
1985  mse_th(itest,ifac) = 1.0/real(nam%ens1_ne-1,kind_real)*sum(1+fld_ref(:,:,:,:,itest)**2)
1986 
1987  ! Write first 10 test vectors
1988  if (itest<=min(ntest,10)) then
1989  ! Write field
1990  write(itestchar,'(i4.4)') itest
1991  filename = trim(nam%prefix)//'_randomize_'//itestchar
1992  do its=1,nam%nts
1993  write(itschar,'(i2.2)') its
1994  do iv=1,nam%nv
1995  call io%fld_write(mpl,nam,geom,filename,trim(nam%varname(iv))//'_rand_'//nechar//'_'//itschar,fld(:,:,iv,its))
1996  end do
1997  end do
1998  end if
1999  end do
2000 
2001  ! Print scores
2002  write(mpl%info,'(a7,a,i4,a,e15.8,a,e15.8)') '','Ensemble size ',nefac(ifac),', MSE (exp. / th.): ', &
2003  & sum(mse(:,ifac))/real(ntest,kind_real),' / ',sum(mse_th(:,ifac))/real(ntest,kind_real)
2004  call flush(mpl%info)
2005 
2006  ! Release memory
2007  call ens%dealloc
2008 end do
2009 
2010 ! Reset namelist variables
2011 nam%ens1_ne = ens1_ne
2012 
2013 end subroutine nicas_test_randomization
2014 
2015 !----------------------------------------------------------------------
2016 ! Subroutine: nicas_test_consistency
2017 ! Purpose: test HDIAG-NICAS consistency with a randomization method
2018 !----------------------------------------------------------------------
2019 subroutine nicas_test_consistency(nicas,mpl,rng,nam,geom,bpar,io,cmat)
2020 
2021 implicit none
2022 
2023 ! Passed variables
2024 class(nicas_type),intent(in) :: nicas ! NICAS data
2025 type(mpl_type),intent(inout) :: mpl ! MPI data
2026 type(rng_type),intent(inout) :: rng ! Random number generator
2027 type(nam_type),intent(inout) :: nam ! Namelist variables
2028 type(geom_type),intent(in) :: geom ! Geometry
2029 type(bpar_type),intent(in) :: bpar ! Block parameters
2030 type(io_type),intent(in) :: io ! I/O
2031 type(cmat_type),intent(in) :: cmat ! C matrix data
2032 
2033 ! Local variables
2034 integer :: ens1_ne,ens1_ne_offset,ens1_nsub,ib,il0
2035 real(kind_real) :: rh_c0_sum,rv_c0_sum,norm
2036 character(len=1024) :: prefix,method
2037 type(cmat_type) :: cmat_test
2038 type(hdiag_type) :: hdiag_test
2039 type(ens_type) :: ens
2040 
2041 ! Randomize ensemble
2042 write(mpl%info,'(a)') '-------------------------------------------------------------------'
2043 write(mpl%info,'(a)') '--- Randomize ensemble'
2044 call flush(mpl%info)
2045 call nicas%randomize(mpl,rng,nam,geom,bpar,ne_rand,ens)
2046 
2047 ! Copy sampling
2048 call execute_command_line('cp -f '//trim(nam%datadir)//'/'//trim(nam%prefix)//'_sampling.nc ' &
2049  & //trim(nam%datadir)//'/'//trim(nam%prefix)//'_consistency-test_sampling.nc')
2050 
2051 ! Save namelist variables
2052 prefix = nam%prefix
2053 method = nam%method
2054 ens1_ne = nam%ens1_ne
2055 ens1_ne_offset = nam%ens1_ne_offset
2056 ens1_nsub = nam%ens1_nsub
2057 
2058 ! Set namelist variables
2059 nam%prefix = trim(nam%prefix)//'_consistency-test'
2060 nam%method = 'cor'
2061 nam%ens1_ne = ne_rand
2062 nam%ens1_ne_offset = 0
2063 nam%ens1_nsub = 1
2064 
2065 ! Call hdiag driver
2066 call hdiag_test%run_hdiag(mpl,rng,nam,geom,bpar,io,ens)
2067 
2068 ! Copy into C matrix
2069 call cmat_test%from_hdiag(mpl,nam,geom,bpar,hdiag_test)
2070 
2071 ! Print scores
2072 write(mpl%info,'(a)') '-------------------------------------------------------------------'
2073 write(mpl%info,'(a)') '--- HDIAG/NICAS consistency results'
2074 do ib=1,bpar%nbe
2075  if (bpar%nicas_block(ib)) then
2076  write(mpl%info,'(a7,a,a)') '','Block: ',trim(bpar%blockname(ib))
2077  do il0=1,geom%nl0
2078  call mpl%f_comm%allreduce(sum(cmat_test%blk(ib)%rh(:,il0)-cmat%blk(ib)%rh(:,il0),geom%mask_c0a(:,il0)), &
2079  & rh_c0_sum,fckit_mpi_sum())
2080  call mpl%f_comm%allreduce(sum(cmat_test%blk(ib)%rv(:,il0)-cmat%blk(ib)%rv(:,il0),geom%mask_c0a(:,il0)),&
2081  & rv_c0_sum,fckit_mpi_sum())
2082  call mpl%f_comm%allreduce(real(count(geom%mask_c0a(:,il0)),kind_real),norm,fckit_mpi_sum())
2083  write(mpl%info,'(a10,a7,i3,a4,a25,f6.1,a)') '','Level: ',nam%levs(il0),' ~> ','horizontal length-scale: ', &
2084  & rh_c0_sum/norm*reqkm,' km'
2085  if (any(abs(cmat%blk(ib)%rv(:,il0))>0.0)) then
2086  write(mpl%info,'(a49,f6.1,a)') 'vertical length-scale: ',rh_c0_sum/norm,' '//trim(mpl%vunitchar)
2087  end if
2088  end do
2089  end if
2090 end do
2091 call flush(mpl%info)
2092 
2093 ! Reset namelist variables
2094 nam%prefix = prefix
2095 nam%method = method
2096 nam%ens1_ne = ens1_ne
2097 nam%ens1_ne_offset = ens1_ne_offset
2098 nam%ens1_nsub = ens1_nsub
2099 
2100 end subroutine nicas_test_consistency
2101 
2102 !----------------------------------------------------------------------
2103 ! Subroutine: nicas_test_optimality
2104 ! Purpose: test HDIAG localization optimality with a randomization method
2105 !----------------------------------------------------------------------
2106 subroutine nicas_test_optimality(nicas,mpl,rng,nam,geom,bpar,io)
2107 
2108 implicit none
2109 
2110 ! Passed variables
2111 class(nicas_type),intent(in) :: nicas ! NICAS data
2112 type(mpl_type),intent(inout) :: mpl ! MPI data
2113 type(rng_type),intent(inout) :: rng ! Random number generator
2114 type(nam_type),intent(inout) :: nam ! Namelist variables
2115 type(geom_type),intent(in) :: geom ! Geometry
2116 type(bpar_type),intent(in) :: bpar ! Block parameters
2117 type(io_type),intent(in) :: io ! I/O
2118 
2119 ! Local variables
2120 integer :: ib,ifac,itest
2121 real(kind_real) :: fld_ref(geom%nc0a,geom%nl0,nam%nv,nam%nts,ntest),fld_save(geom%nc0a,geom%nl0,nam%nv,nam%nts,ntest)
2122 real(kind_real) :: fld(geom%nc0a,geom%nl0,nam%nv,nam%nts),fac(nfac),mse(ntest,nfac)
2123 character(len=1024) :: prefix,method
2124 type(cmat_type) :: cmat_save,cmat_test
2125 type(hdiag_type) :: hdiag_save
2126 type(ens_type) :: ens
2127 type(nicas_type) :: nicas_test
2128 
2129 ! Define test vectors
2130 write(mpl%info,'(a4,a)') '','Define test vectors'
2131 call flush(mpl%info)
2132 call define_test_vectors(mpl,rng,nam,geom,ntest,fld_save)
2133 
2134 ! Apply NICAS to test vectors
2135 write(mpl%info,'(a4,a)') '','Apply NICAS to test vectors'
2136 call flush(mpl%info)
2137 fld_ref = fld_save
2138 do itest=1,ntest
2139  call nicas%apply_from_sqrt(mpl,nam,geom,bpar,fld_ref(:,:,:,:,itest))
2140 end do
2141 
2142 ! Randomize ensemble
2143 write(mpl%info,'(a)') '-------------------------------------------------------------------'
2144 write(mpl%info,'(a)') '--- Randomize ensemble'
2145 call flush(mpl%info)
2146 call nicas%randomize(mpl,rng,nam,geom,bpar,nam%ens1_ne,ens)
2147 
2148 ! Copy sampling
2149 call execute_command_line('cp -f '//trim(nam%datadir)//'/'//trim(nam%prefix)//'_sampling.nc ' &
2150  & //trim(nam%datadir)//'/'//trim(nam%prefix)//'_optimality-test_sampling.nc')
2151 
2152 ! Save namelist variables
2153 prefix = nam%prefix
2154 method = nam%method
2155 
2156 ! Set namelist variables
2157 nam%prefix = trim(nam%prefix)//'_optimality-test'
2158 nam%method = 'loc_norm'
2159 
2160 ! Allocation
2161 call nicas_test%alloc(mpl,nam,bpar,'nicas_test')
2162 
2163 ! Call hdiag driver
2164 call hdiag_save%run_hdiag(mpl,rng,nam,geom,bpar,io,ens)
2165 
2166 ! Copy into C matrix
2167 call cmat_save%from_hdiag(mpl,nam,geom,bpar,hdiag_save)
2168 
2169 ! Copy cmat
2170 cmat_test = cmat_save%copy(nam,geom,bpar)
2171 
2172 do ifac=1,nfac
2173  ! Multiplication factor
2174  fac(ifac) = 2.0*real(ifac,kind_real)/real(nfac,kind_real)
2175 
2176  write(mpl%info,'(a)') '-------------------------------------------------------------------'
2177  write(mpl%info,'(a,f4.2,a)') '--- Apply a multiplicative factor ',fac(ifac),' to length-scales'
2178  call flush(mpl%info)
2179 
2180  do ib=1,bpar%nbe
2181  if (bpar%nicas_block(ib)) then
2182  ! Length-scales multiplication
2183  cmat_test%blk(ib)%rh = fac(ifac)*cmat_save%blk(ib)%rh
2184  cmat_test%blk(ib)%rv = fac(ifac)*cmat_save%blk(ib)%rv
2185  if (trim(nam%strategy)=='specific_multivariate') then
2186  cmat_test%blk(ib)%rhs = fac(ifac)*cmat_save%blk(ib)%rhs
2187  cmat_test%blk(ib)%rvs = fac(ifac)*cmat_save%blk(ib)%rvs
2188  end if
2189 
2190  ! Compute NICAS parameters
2191  call nicas_test%blk(ib)%compute_parameters(mpl,rng,nam,geom,cmat_test%blk(ib))
2192  end if
2193 
2194  if (bpar%B_block(ib)) then
2195  ! Copy weights
2196  nicas_test%blk(ib)%wgt = cmat_test%blk(ib)%wgt
2197  if (bpar%nicas_block(ib)) then
2198  allocate(nicas_test%blk(ib)%coef_ens(geom%nc0a,geom%nl0))
2199  nicas_test%blk(ib)%coef_ens = cmat_test%blk(ib)%coef_ens
2200  end if
2201  end if
2202  end do
2203 
2204  do itest=1,ntest
2205  ! Test NICAS
2206  fld = fld_save(:,:,:,:,itest)
2207  call nicas_test%apply_bens(mpl,nam,geom,bpar,ens,fld)
2208 
2209  ! RMSE
2210  mse(itest,ifac) = sum((fld-fld_ref(:,:,:,:,itest))**2)
2211  end do
2212 
2213  ! Print scores
2214  write(mpl%info,'(a)') '-------------------------------------------------------------------'
2215  write(mpl%info,'(a,f4.2,a,e15.8)') '--- Optimality results for a factor ',fac(ifac),', MSE: ', &
2216  & sum(mse(:,ifac))/real(ntest,kind_real)
2217  call flush(mpl%info)
2218 end do
2219 
2220 ! Print scores summary
2221 write(mpl%info,'(a)') '-------------------------------------------------------------------'
2222 write(mpl%info,'(a)') '--- Optimality results summary'
2223 do ifac=1,nfac
2224  write(mpl%info,'(a7,a,f4.2,a,e15.8)') '','Factor ',fac(ifac),', MSE: ',sum(mse(:,ifac))/real(ntest,kind_real)
2225 end do
2226 call flush(mpl%info)
2227 
2228 ! Reset namelist variables
2229 nam%prefix = prefix
2230 nam%method = method
2231 
2232 end subroutine nicas_test_optimality
2233 
2234 !----------------------------------------------------------------------
2235 ! Subroutine: define_test_vectors
2236 ! Purpose: define test vectors
2237 !----------------------------------------------------------------------
2238 subroutine define_test_vectors(mpl,rng,nam,geom,ntest,fld)
2239 
2240 ! Passed variables
2241 type(mpl_type),intent(in) :: mpl ! MPI data
2242 type(rng_type),intent(inout) :: rng ! Random number generator
2243 type(nam_type),intent(in) :: nam ! Namelist
2244 type(geom_type),intent(in) :: geom ! Geometry
2245 integer,intent(in) :: ntest ! Number of vectors
2246 real(kind_real),intent(out) :: fld(geom%nc0a,geom%nl0,nam%nv,nam%nts,ntest) ! Field
2247 
2248 ! Local variables
2249 integer :: ic0dir(ntest),il0dir(ntest),ivdir(ntest),itsdir(ntest)
2250 integer :: itest,ic0,iproc,ic0a
2251 
2252 ! Define random dirac locations
2253 if (mpl%main) then
2254  call rng%rand_integer(1,geom%nc0,ic0dir)
2255  call rng%rand_integer(1,geom%nl0,il0dir)
2256 end if
2257 call mpl%f_comm%broadcast(ic0dir,mpl%ioproc-1)
2258 call mpl%f_comm%broadcast(il0dir,mpl%ioproc-1)
2259 ivdir = 1
2260 itsdir = 1
2261 
2262 ! Define test vectors
2263 do itest=1,ntest
2264  fld(:,:,:,:,itest) = 0.0
2265  ic0 = ic0dir(itest)
2266  iproc = geom%c0_to_proc(ic0)
2267  if (iproc==mpl%myproc) then
2268  ic0a = geom%c0_to_c0a(ic0)
2269  fld(ic0a,il0dir(itest),ivdir(itest),itsdir(itest),itest) = 1.0
2270  end if
2271 end do
2272 
2273 end subroutine define_test_vectors
2274 
2275 end module type_nicas
type_nicas::nicas_test_dirac
subroutine nicas_test_dirac(nicas, mpl, nam, geom, bpar, io, ens)
Definition: type_nicas.F90:1854
type_nicas::nicas_alloc_cv
subroutine nicas_alloc_cv(nicas, bpar, cv, getsizeonly)
Definition: type_nicas.F90:674
type_nicas::nicas_type
Definition: type_nicas.F90:39
type_nicas::nicas_randomize
subroutine nicas_randomize(nicas, mpl, rng, nam, geom, bpar, ne, ens)
Definition: type_nicas.F90:1472
type_rng
Definition: type_rng.F90:8
type_nicas::pos_def_test
logical, parameter pos_def_test
Definition: type_nicas.F90:36
type_nicas::nfac
integer, parameter nfac
Definition: type_nicas.F90:34
type_nicas::nicas_read
subroutine nicas_read(nicas, mpl, nam, geom, bpar)
Definition: type_nicas.F90:145
tools_const::rad2deg
real(kind_real), parameter, public rad2deg
Definition: tools_const.F90:17
tools_const::msvali
integer, parameter, public msvali
Definition: tools_const.F90:23
type_nicas::nicas_test_consistency
subroutine nicas_test_consistency(nicas, mpl, rng, nam, geom, bpar, io, cmat)
Definition: type_nicas.F90:2020
type_nicas::nicas_run_nicas_tests
subroutine nicas_run_nicas_tests(nicas, mpl, rng, nam, geom, bpar, io, cmat, ens)
Definition: type_nicas.F90:544
type_cv
Definition: type_cv.F90:8
tools_nc
Definition: tools_nc.F90:8
type_nicas
Definition: type_nicas.F90:8
type_bpar
Definition: type_bpar.F90:8
type_nicas::ntest
integer, parameter ntest
Definition: type_nicas.F90:35
type_nicas_blk
Definition: type_nicas_blk.F90:8
type_nicas_blk::nicas_blk_type
Definition: type_nicas_blk.F90:58
type_cv::cv_type
Definition: type_cv.F90:16
type_cmat
Definition: type_cmat.F90:8
type_nicas::nicas_test_optimality
subroutine nicas_test_optimality(nicas, mpl, rng, nam, geom, bpar, io)
Definition: type_nicas.F90:2107
tools_kinds
Definition: tools_kinds.F90:8
tools_func::cholesky
subroutine, public cholesky(mpl, n, a, u)
Definition: tools_func.F90:745
type_nicas::nicas_random_cv
subroutine nicas_random_cv(nicas, rng, bpar, cv)
Definition: type_nicas.F90:726
type_cmat::cmat_type
Definition: type_cmat.F90:35
type_nam::nam_type
Definition: type_nam.F90:29
tools_const::msvalr
real(kind_real), parameter, public msvalr
Definition: tools_const.F90:24
type_nicas::nicas_apply_from_sqrt
subroutine nicas_apply_from_sqrt(nicas, mpl, nam, geom, bpar, fld)
Definition: type_nicas.F90:1008
type_nicas::ne_rand
integer, parameter ne_rand
Definition: type_nicas.F90:33
type_hdiag::hdiag_type
Definition: type_hdiag.F90:27
type_nam
Definition: type_nam.F90:8
type_nicas::nicas_write_mpi_summary
subroutine nicas_write_mpi_summary(nicas, mpl, nam, geom, bpar)
Definition: type_nicas.F90:395
type_nicas::nicas_apply_sqrt_ad
subroutine nicas_apply_sqrt_ad(nicas, mpl, nam, geom, bpar, fld, cv)
Definition: type_nicas.F90:1251
tools_const
Definition: tools_const.F90:8
type_nicas::nicas_dealloc
subroutine nicas_dealloc(nicas, nam, geom, bpar)
Definition: type_nicas.F90:117
tools_func
Definition: tools_func.F90:8
type_mpl::mpl_type
Definition: type_mpl.F90:25
type_rng::rng_type
Definition: type_rng.F90:29
type_linop::linop_type
Definition: type_linop.F90:32
tools_missing::msr
Definition: tools_missing.F90:24
type_mpl
Definition: type_mpl.F90:8
type_nicas::nicas_write
subroutine nicas_write(nicas, mpl, nam, geom, bpar)
Definition: type_nicas.F90:288
type_io
Definition: type_io.F90:8
type_nicas::nicas_apply
subroutine nicas_apply(nicas, mpl, nam, geom, bpar, fld)
Definition: type_nicas.F90:753
type_com
Definition: type_com.F90:8
tools_missing::msi
Definition: tools_missing.F90:15
type_nicas::nicas_apply_bens_noloc
subroutine nicas_apply_bens_noloc(nicas, mpl, nam, geom, ens, fld)
Definition: type_nicas.F90:1602
max
#define max(a, b)
Definition: mosaic_util.h:33
type_linop
Definition: type_linop.F90:8
type_nicas::nicas_test_adjoint
subroutine nicas_test_adjoint(nicas, mpl, rng, nam, geom, bpar, ens)
Definition: type_nicas.F90:1657
type_nicas::nicas_alloc
subroutine nicas_alloc(nicas, mpl, nam, bpar, prefix)
Definition: type_nicas.F90:78
type_geom
Definition: type_geom.F90:8
type_nicas::nicas_apply_sqrt
subroutine nicas_apply_sqrt(nicas, mpl, nam, geom, bpar, cv, fld)
Definition: type_nicas.F90:1050
type_nicas::nicas_test_sqrt
subroutine nicas_test_sqrt(nicas, mpl, rng, nam, geom, bpar, io, cmat, ens)
Definition: type_nicas.F90:1739
type_geom::geom_type
Definition: type_geom.F90:33
type_ens::ens_type
Definition: type_ens.F90:19
type_nicas::nicas_test_randomization
subroutine nicas_test_randomization(nicas, mpl, rng, nam, geom, bpar, io)
Definition: type_nicas.F90:1917
tools_missing
Definition: tools_missing.F90:8
fv3jedi_kinds_mod::kind_real
integer, parameter, public kind_real
Definition: fv3jedi_kinds_mod.f90:13
min
#define min(a, b)
Definition: mosaic_util.h:32
type_io::io_type
Definition: type_io.F90:27
type_nicas::define_test_vectors
subroutine define_test_vectors(mpl, rng, nam, geom, ntest, fld)
Definition: type_nicas.F90:2239
type_nicas::nicas_run_nicas
subroutine nicas_run_nicas(nicas, mpl, rng, nam, geom, bpar, cmat)
Definition: type_nicas.F90:478
type_com::com_type
Definition: type_com.F90:21
tools_nc::ncfloat
integer, parameter, public ncfloat
Definition: tools_nc.F90:27
type_ens
Definition: type_ens.F90:8
type_bpar::bpar_type
Definition: type_bpar.F90:17
type_hdiag
Definition: type_hdiag.F90:8
type_nicas::nicas_apply_bens
subroutine nicas_apply_bens(nicas, mpl, nam, geom, bpar, ens, fld)
Definition: type_nicas.F90:1543
tools_missing::isnotmsr
Definition: tools_missing.F90:49
crtm_parameters::pi
real(fp), parameter, public pi
Definition: CRTM_Parameters.f90:86
tools_const::reqkm
real(kind_real), parameter, public reqkm
Definition: tools_const.F90:19