ROMS
Loading...
Searching...
No Matches
correlation.h
Go to the documentation of this file.
1 MODULE roms_kernel_mod
2!
3!git $Id$
4!================================================== Hernan G. Arango ===
5! Copyright (c) 2002-2025 The ROMS Group !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md !
8!=======================================================================
9! !
10! ROMS 4DVAR Background-error Correlation Model !
11! !
12! This driver is used to build and test the 4DVAR background-error !
13! correlation model using a generalized difussion operator: !
14! !
15! B = S C S !
16! !
17! C = C^(1/2) C^(T/2) !
18! !
19! C^(1/2) = G L^(1/2) W^(-1/2) !
20! !
21! C^(T/2) = W^(T/2) L^(T/2) G !
22! !
23! where !
24! !
25! B : background-error covariance !
26! C : background-error correlation !
27! G : normalization coefficient matrix !
28! L : self-adjoint diffusion filter !
29! S : background-error standard deviation !
30! W : Grid cell area or volume diagonal matrix !
31! !
32! The routines in this driver control the initialization, time- !
33! stepping, and finalization of ROMS model following ESMF !
34! conventions: !
35! !
36! ROMS_initialize !
37! ROMS_run !
38! ROMS_finalize !
39! !
40! Reference !
41! !
42! Weaver, A. and P. Courtier, 2001: Correlation modelling on !
43! the sphere using a generalized diffusion equation, Q. J. !
44! R. Meteorol. Soc., 127, 1815-1845. !
45! !
46!=======================================================================
47!
48 USE mod_param
49 USE mod_parallel
50 USE mod_arrays
51 USE mod_fourdvar
52 USE mod_iounits
53 USE mod_ncparam
54 USE mod_scalars
55 USE mod_stepping
56!
57#ifdef BALANCE_OPERATOR
58 USE ad_balance_mod, ONLY : ad_balance
59#endif
61 USE ad_def_his_mod, ONLY : ad_def_his
63 USE ad_wrt_his_mod, ONLY : ad_wrt_his
64 USE analytical_mod, ONLY : ana_perturb
66 USE def_norm_mod, ONLY : def_norm
67 USE get_state_mod, ONLY : get_state
68 USE ini_adjust_mod, ONLY : load_adtotl
69 USE ini_adjust_mod, ONLY : load_tltoad
70 USE inp_par_mod, ONLY : inp_par
73 USE strings_mod, ONLY : founderror
74#ifdef BALANCE_OPERATOR
75 USE tl_balance_mod, ONLY : tl_balance
76#endif
79 USE strings_mod, ONLY : founderror
80 USE wrt_rst_mod, ONLY : wrt_rst
81#if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
83#endif
84!
85 implicit none
86!
87 PUBLIC :: roms_initialize
88 PUBLIC :: roms_run
89 PUBLIC :: roms_finalize
90!
91 CONTAINS
92!
93 SUBROUTINE roms_initialize (first, mpiCOMM)
94!
95!=======================================================================
96! !
97! This routine allocates and initializes ROMS state variables !
98! and internal and external parameters. !
99! !
100!=======================================================================
101!
102! Imported variable declarations.
103!
104 logical, intent(inout) :: first
105!
106 integer, intent(in), optional :: mpiCOMM
107!
108! Local variable declarations.
109!
110 logical :: allocate_vars = .true.
111!
112#ifdef DISTRIBUTE
113 integer :: MyError, MySize
114#endif
115 integer :: STDrec, Tindex
116 integer :: chunk_size, ng, thread, tile
117#ifdef _OPENMP
118 integer :: my_threadnum
119#endif
120!
121 character (len=*), parameter :: MyFile = &
122 & __FILE__//", ROMS_initialize"
123
124#ifdef DISTRIBUTE
125!
126!-----------------------------------------------------------------------
127! Set distribute-memory (mpi) world communicator.
128!-----------------------------------------------------------------------
129!
130 IF (PRESENT(mpicomm)) THEN
131 ocn_comm_world=mpicomm
132 ELSE
133 ocn_comm_world=mpi_comm_world
134 END IF
135 CALL mpi_comm_rank (ocn_comm_world, myrank, myerror)
136 CALL mpi_comm_size (ocn_comm_world, mysize, myerror)
137#endif
138!
139!-----------------------------------------------------------------------
140! On first pass, initialize model parameters a variables for all
141! nested/composed grids. Notice that the logical switch "first"
142! is used to allow multiple calls to this routine during ensemble
143! configurations.
144!-----------------------------------------------------------------------
145!
146 IF (first) THEN
147 first=.false.
148!
149! Initialize parallel control switches. These scalars switches are
150! independent from standard input parameters.
151!
153!
154! Set the ROMS standard output unit to write verbose execution info.
155! Notice that the default standard out unit in Fortran is 6.
156!
157! In some applications like coupling or disjointed mpi-communications,
158! it is advantageous to write standard output to a specific filename
159! instead of the default Fortran standard output unit 6. If that is
160! the case, it opens such formatted file for writing.
161!
162 IF (set_stdoutunit) THEN
164 set_stdoutunit=.false.
165 END IF
166!
167! Read in model tunable parameters from standard input. Allocate and
168! initialize variables in several modules after the number of nested
169! grids and dimension parameters are known.
170!
171 CALL inp_par (inlm)
172 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
173!
174! Set domain decomposition tile partition range. This range is
175! computed only once since the "first_tile" and "last_tile" values
176! are private for each parallel thread/node.
177!
178#if defined _OPENMP
179 mythread=my_threadnum()
180#elif defined DISTRIBUTE
182#else
183 mythread=0
184#endif
185 DO ng=1,ngrids
186 chunk_size=(ntilex(ng)*ntilee(ng)+numthreads-1)/numthreads
187 first_tile(ng)=mythread*chunk_size
188 last_tile(ng)=first_tile(ng)+chunk_size-1
189 END DO
190!
191! Initialize internal wall clocks. Notice that the timings does not
192! includes processing standard input because several parameters are
193! needed to allocate clock variables.
194!
195 IF (master) THEN
196 WRITE (stdout,10)
197 10 FORMAT (/,' Process Information:',/)
198 END IF
199!
200 DO ng=1,ngrids
201 DO thread=thread_range
202 CALL wclock_on (ng, inlm, 0, __line__, myfile)
203 END DO
204 END DO
205!
206! Allocate and initialize modules variables.
207!
208 CALL roms_allocate_arrays (allocate_vars)
210 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
211
212 END IF
213!
214!-----------------------------------------------------------------------
215! Initialize metrics over all nested grids, if applicable.
216!-----------------------------------------------------------------------
217!
218 CALL initial
219 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
220!
221! Adjust "time" variable since we are not time-stepping.
222!
223 DO ng=1,ngrids
224 time(ng)=time(ng)+dt(ng)
225 END DO
226!
227! Initialize run or ensemble counter.
228!
229 nrun=1
230!
231!-----------------------------------------------------------------------
232! Read in standard deviation factors for error covariance.
233!-----------------------------------------------------------------------
234!
235! Initial conditions standard deviation. They are loaded in Tindex=1
236! of the e_var(...,Tindex) state variables.
237!
238 stdrec=1
239 tindex=1
240 DO ng=1,ngrids
241 IF (ldefnrm(1,ng).or.lwrtnrm(1,ng)) THEN
242 CALL get_state (ng, 10, 10, std(1,ng), stdrec, tindex)
243 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
244 END IF
245 END DO
246!
247! Model error standard deviation. They are loaded in Tindex=2
248! of the e_var(...,Tindex) state variables.
249!
250 stdrec=1
251 tindex=2
252 DO ng=1,ngrids
253 IF ((ldefnrm(2,ng).or.lwrtnrm(2,ng)).and.(nsa.eq.2)) THEN
254 CALL get_state (ng, 11, 11, std(2,ng), stdrec, tindex)
255 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
256 END IF
257 END DO
258
259#ifdef ADJUST_BOUNDARY
260!
261! Open boundary conditions standard deviation.
262!
263 stdrec=1
264 tindex=1
265 DO ng=1,ngrids
266 IF (ldefnrm(3,ng).or.lwrtnrm(3,ng)) THEN
267 CALL get_state (ng, 12, 12, std(3,ng), stdrec, tindex)
268 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
269 END IF
270 END DO
271#endif
272#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
273!
274! Surface forcing standard deviation.
275!
276 stdrec=1
277 tindex=1
278 DO ng=1,ngrids
279 IF (ldefnrm(4,ng).or.lwrtnrm(4,ng)) THEN
280 CALL get_state (ng, 13, 13, std(4,ng), stdrec, tindex)
281 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
282 END IF
283 END DO
284#endif
285!
286!-----------------------------------------------------------------------
287! Compute or read in error covariance normalization factors.
288!-----------------------------------------------------------------------
289!
290! If computing, write out factors to NetCDF. This is an expensive
291! computation and needs to be computed once for a particular
292! application grid.
293!
294 DO ng=1,ngrids
295 IF (any(lwrtnrm(:,ng))) THEN
296 IF (ldefnrm(1,ng).or.lwrtnrm(1,ng)) THEN
297 CALL def_norm (ng, inlm, 1)
298 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
299 END IF
300
301 IF ((ldefnrm(2,ng).or.lwrtnrm(2,ng)).and.(nsa.eq.2)) THEN
302 CALL def_norm (ng, inlm, 2)
303 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
304 END IF
305#ifdef ADJUST_BOUNDARY
306 IF (ldefnrm(3,ng).or.lwrtnrm(3,ng)) THEN
307 CALL def_norm (ng, inlm, 3)
308 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
309 END IF
310#endif
311#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
312 IF (ldefnrm(4,ng).or.lwrtnrm(4,ng)) THEN
313 CALL def_norm (ng, inlm, 4)
314 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
315 END IF
316#endif
317 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
318 DO tile=first_tile(ng),last_tile(ng),+1
319 CALL normalization (ng, tile, 2)
320 END DO
321 ldefnrm(1:4,ng)=.false.
322 lwrtnrm(1:4,ng)=.false.
323 END IF
324 END DO
325!
326 RETURN
327 END SUBROUTINE roms_initialize
328!
329 SUBROUTINE roms_run (RunInterval)
330!
331!=======================================================================
332! !
333! This routine computes background-error correlations. !
334! !
335!=======================================================================
336!
337! Imported variable declarations.
338!
339 real(dp), intent(in) :: RunInterval ! seconds
340!
341! Local variable declarations.
342!
343 logical :: Lweak, add
344 integer :: i, ng, tile
345#ifdef BALANCE_OPERATOR
346 integer :: Lbck = 1
347#endif
348!
349 character (len=*), parameter :: MyFile = &
350 & __FILE__//", ROMS_run"
351!
352!-----------------------------------------------------------------------
353! Test correlation model.
354!-----------------------------------------------------------------------
355
356#ifdef BALANCE_OPERATOR
357!
358! Read background state, use initial conditions.
359!
360 DO ng=1,ngrids
361 CALL get_state (ng, inlm, 9, ini(ng), lbck, lbck)
362 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
363 END DO
364
365# ifdef ZETA_ELLIPTIC
366!
367! Compute the reference zeta and biconjugate gradient arrays
368! required for the balance of free surface.
369!
370 IF (balance(isfsur)) THEN
371 DO ng=1,ngrids
372 DO tile=first_tile(ng),last_tile(ng),+1
373 CALL balance_ref (ng, tile, lbck)
374 CALL biconj (ng, tile, inlm, lbck)
375 END DO
376 wrtzetaref(ng)=.true.
377 END DO
378 END IF
379# endif
380#endif
381!
382! Initialize adjoint model state with a delta function at specified
383! point. Use USER parameters from standard input to perturb solution
384! in routine "ana_perturb". Then, convolve solution with the adjoint
385! diffusion operator.
386!
387 admodel=.true.
388 lweak=.false.
389
390 DO ng=1,ngrids
391 lnew(ng)=1
392 DO tile=first_tile(ng),last_tile(ng),+1
393 CALL ana_perturb (ng, tile, iadm)
394#ifdef BALANCE_OPERATOR
395 CALL ad_balance (ng, tile, lbck, lnew(ng))
396 CALL ad_variability (ng, tile, lnew(ng), lweak)
397#endif
398 CALL ad_convolution (ng, tile, lnew(ng), lweak, 2)
399 END DO
400
401 admodel=.false.
402!
403! Initialize tangent linear model with convolved adjoint solution.
404! Then, apply tangent linear convolution.
405!
406 add=.false.
407 DO tile=first_tile(ng),last_tile(ng),+1
408 CALL load_adtotl (ng, tile, lnew(ng), lnew(ng), add)
409 CALL tl_convolution (ng, tile, lnew(ng), lweak, 2)
410#ifdef BALANCE_OPERATOR
411 CALL tl_variability (ng, tile, lnew(ng), lweak)
412 CALL tl_balance (ng, tile, lbck, lnew(ng))
413#endif
414 CALL load_tltoad (ng, tile, lnew(ng), lnew(ng), add)
415 END DO
416 END DO
417!
418! Write out background error correlation in adjoint history NetCDF
419! file.
420!
421 DO ng=1,ngrids
422 kstp(ng)=lnew(ng)
423#ifdef SOLVE3D
424 nstp(ng)=lnew(ng)
425#endif
426 ldefadj(ng)=.true.
427 lwrtadj(ng)=.true.
428 lwrtstate2d(ng)=.true.
429 CALL ad_def_his (ng, ldefadj(ng))
430 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
431#if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
432 ladjusted(ng)=.true.
433#endif
434#ifdef DISTRIBUTE
435 CALL ad_wrt_his (ng, myrank)
436#else
437 CALL ad_wrt_his (ng, -1)
438#endif
439 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
440#if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
441 ladjusted(ng)=.false.
442#endif
443 END DO
444!
445 RETURN
446 END SUBROUTINE roms_run
447!
448 SUBROUTINE roms_finalize
449!
450!=======================================================================
451! !
452! This routine terminates ROMS driver execution. !
453! !
454!=======================================================================
455!
456! Local variable declarations.
457!
458 integer :: Fcount, ng, thread
459!
460 character (len=*), parameter :: MyFile = &
461 & __FILE__//", ROMS_finalize"
462!
463!-----------------------------------------------------------------------
464! If blowing-up, save latest model state into RESTART NetCDF file.
465!-----------------------------------------------------------------------
466!
467! If cycling restart records, write solution into record 3.
468!
469 IF (exit_flag.eq.1) THEN
470 DO ng=1,ngrids
471 IF (lwrtrst(ng)) THEN
472 IF (master) WRITE (stdout,10)
473 10 FORMAT (/,' Blowing-up: Saving latest model state into ', &
474 & ' RESTART file',/)
475 fcount=rst(ng)%load
476 IF (lcyclerst(ng).and.(rst(ng)%Nrec(fcount).ge.2)) THEN
477 rst(ng)%Rindex=2
478 lcyclerst(ng)=.false.
479 END IF
482#ifdef DISTRIBUTE
483 CALL wrt_rst (ng, myrank)
484#else
485 CALL wrt_rst (ng, -1)
486#endif
487 END IF
488 END DO
489 END IF
490!
491!-----------------------------------------------------------------------
492! Stop model and time profiling clocks, report memory requirements, and
493! close output NetCDF files.
494!-----------------------------------------------------------------------
495!
496! Stop time clocks.
497!
498 IF (master) THEN
499 WRITE (stdout,20)
500 20 FORMAT (/,'Elapsed wall CPU time for each process (seconds):',/)
501 END IF
502!
503 DO ng=1,ngrids
504 DO thread=thread_range
505 CALL wclock_off (ng, inlm, 0, __line__, myfile)
506 END DO
507 END DO
508!
509! Report dynamic memory and automatic memory requirements.
510!
511 CALL memory
512!
513! Close IO files.
514!
515 DO ng=1,ngrids
516 CALL close_inp (ng, inlm)
517 END DO
518 CALL close_out
519!
520 RETURN
521 END SUBROUTINE roms_finalize
522
523 END MODULE roms_kernel_mod
subroutine initial
Definition initial.F:3
subroutine memory
Definition memory.F:3
subroutine, public ad_balance(ng, tile, lbck, linp)
Definition ad_balance.F:59
subroutine, public ad_convolution(ng, tile, linp, lweak, ifac)
subroutine, public ad_def_his(ng, ldef)
Definition ad_def_his.F:51
subroutine, public ad_variability(ng, tile, linp, lweak)
subroutine, public ad_wrt_his(ng, tile)
Definition ad_wrt_his.F:66
subroutine ana_perturb(ng, tile, model)
Definition ana_perturb.h:3
subroutine, public close_out
Definition close_io.F:175
subroutine, public close_inp(ng, model)
Definition close_io.F:92
subroutine, public def_norm(ng, model, ifile)
Definition def_norm.F:53
subroutine, public get_state(ng, model, msg, s, inirec, tindex)
Definition get_state.F:90
subroutine, public load_adtotl(ng, tile, linp, lout, add)
Definition ini_adjust.F:814
subroutine, public load_tltoad(ng, tile, linp, lout, add)
subroutine, public inp_par(model)
Definition inp_par.F:56
subroutine, public roms_initialize_arrays
Definition mod_arrays.F:351
subroutine, public roms_allocate_arrays(allocate_vars)
Definition mod_arrays.F:114
logical, dimension(:), allocatable wrtzetaref
type(t_io), dimension(:,:), allocatable std
type(t_io), dimension(:), allocatable rst
type(t_io), dimension(:), allocatable ini
integer stdout
integer isfsur
subroutine, public initialize_parallel
integer numthreads
integer, dimension(:), allocatable first_tile
integer mythread
logical master
integer, dimension(:), allocatable last_tile
integer ocn_comm_world
integer, parameter inlm
Definition mod_param.F:662
integer, dimension(:), allocatable ntilex
Definition mod_param.F:685
integer, parameter iadm
Definition mod_param.F:665
integer ngrids
Definition mod_param.F:113
integer, dimension(:), allocatable ntilee
Definition mod_param.F:686
integer nsa
Definition mod_param.F:612
logical, dimension(:,:), allocatable lwrtnrm
logical, dimension(:), allocatable ladjusted
real(dp), dimension(:), allocatable dt
integer blowup
logical, dimension(:), allocatable balance
logical admodel
logical, dimension(:), allocatable ldefadj
logical, dimension(:), allocatable lwrtadj
logical, dimension(:,:), allocatable ldefnrm
integer exit_flag
logical, dimension(:), allocatable lwrtstate2d
logical, dimension(:), allocatable lwrtrst
real(dp), dimension(:), allocatable time
integer nrun
integer noerror
logical, dimension(:), allocatable lcyclerst
integer, dimension(:), allocatable kstp
integer, dimension(:), allocatable lnew
integer, dimension(:), allocatable nstp
subroutine, public normalization(ng, tile, ifac)
subroutine, public roms_finalize
Definition ad_roms.h:283
subroutine, public roms_run(runinterval)
Definition ad_roms.h:239
subroutine, public roms_initialize(first, mpicomm)
Definition ad_roms.h:52
integer function, public stdout_unit(mymaster)
Definition stdout_mod.F:48
logical, save set_stdoutunit
Definition stdout_mod.F:41
logical function, public founderror(flag, noerr, line, routine)
Definition strings.F:52
subroutine, public tl_balance(ng, tile, lbck, linp)
Definition tl_balance.F:59
subroutine, public tl_convolution(ng, tile, linp, lweak, ifac)
subroutine, public tl_variability(ng, tile, linp, lweak)
subroutine, public wrt_rst(ng, tile)
Definition wrt_rst.F:63
subroutine, public biconj(ng, tile, model, lbck)
subroutine, public balance_ref(ng, tile, lbck)
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3