SW_VAR_DENS_MODEL  0.9
Dept-averagedgas-particlemodel
solver_2d.f90
Go to the documentation of this file.
1 !********************************************************************************
3 !
6 !
10 !
11 !********************************************************************************
12 MODULE solver_2d
13 
14  ! external variables
15 
16  USE constitutive_2d, ONLY : implicit_flag
17 
18  USE geometry_2d, ONLY : comp_cells_x,comp_cells_y,comp_cells_xy
19  USE geometry_2d, ONLY : comp_interfaces_x,comp_interfaces_y
20 
21  USE geometry_2d, ONLY : b_cent , b_prime_x , b_prime_y
22  USE geometry_2d, ONLY : b_interfacel , b_interfacer , b_interfacet , b_interfaceb
23  USE geometry_2d, ONLY : grav_surf
24  USE geometry_2d, ONLY : source_cell
25 
26  USE parameters_2d, ONLY : n_eqns , n_vars , n_nh , n_solid
27  USE parameters_2d, ONLY : n_rk
28  USE parameters_2d, ONLY : verbose_level
29  USE parameters_2d, ONLY : radial_source_flag
30 
31  USE parameters_2d, ONLY : bcw , bce , bcs , bcn
32 
33  ! external procedures
34  USE geometry_2d, ONLY : limit
35 
36  IMPLICIT none
37 
39  REAL*8 :: t
40 
42  REAL*8, ALLOCATABLE :: q(:,:,:)
44  REAL*8, ALLOCATABLE :: q0(:,:,:)
46  REAL*8, ALLOCATABLE :: q_fv(:,:,:)
47 
49  REAL*8, ALLOCATABLE :: q_interfacel(:,:,:)
51  REAL*8, ALLOCATABLE :: q_interfacer(:,:,:)
53  REAL*8, ALLOCATABLE :: q_interfaceb(:,:,:)
55  REAL*8, ALLOCATABLE :: q_interfacet(:,:,:)
56 
58  REAL*8, ALLOCATABLE :: qp_interfacel(:,:,:)
60  REAL*8, ALLOCATABLE :: qp_interfacer(:,:,:)
62  REAL*8, ALLOCATABLE :: qp_interfaceb(:,:,:)
64  REAL*8, ALLOCATABLE :: qp_interfacet(:,:,:)
65 
67  REAL*8, ALLOCATABLE :: q_cellnw(:,:,:)
69  REAL*8, ALLOCATABLE :: q_cellne(:,:,:)
71  REAL*8, ALLOCATABLE :: q_cellsw(:,:,:)
73  REAL*8, ALLOCATABLE :: q_cellse(:,:,:)
74 
76  REAL*8, ALLOCATABLE :: qp_cellnw(:,:,:)
78  REAL*8, ALLOCATABLE :: qp_cellne(:,:,:)
80  REAL*8, ALLOCATABLE :: qp_cellsw(:,:,:)
82  REAL*8, ALLOCATABLE :: qp_cellse(:,:,:)
83 
84 
86  REAL*8, ALLOCATABLE :: q1max(:,:)
87 
89  REAL*8, ALLOCATABLE :: a_interface_x_max(:,:,:)
91  REAL*8, ALLOCATABLE :: a_interface_y_max(:,:,:)
92 
93 
95  REAL*8, ALLOCATABLE :: a_interface_xneg(:,:,:)
97  REAL*8, ALLOCATABLE :: a_interface_xpos(:,:,:)
99  REAL*8, ALLOCATABLE :: a_interface_yneg(:,:,:)
101  REAL*8, ALLOCATABLE :: a_interface_ypos(:,:,:)
103  REAL*8, ALLOCATABLE :: h_interface_x(:,:,:)
105  REAL*8, ALLOCATABLE :: h_interface_y(:,:,:)
107  REAL*8, ALLOCATABLE :: qp(:,:,:)
108 
110  REAL*8, ALLOCATABLE :: source_xy(:,:)
111 
112  LOGICAL, ALLOCATABLE :: solve_mask(:,:)
113  LOGICAL, ALLOCATABLE :: solve_mask_x(:,:)
114  LOGICAL, ALLOCATABLE :: solve_mask_y(:,:)
115 
116  INTEGER :: solve_cells
117  INTEGER :: solve_interfaces_x
118  INTEGER :: solve_interfaces_y
119 
121  REAL*8 :: dt
122 
123  LOGICAL, ALLOCATABLE :: mask22(:,:) , mask21(:,:) , mask11(:,:) , mask12(:,:)
124 
125  INTEGER :: i_rk
126 
128  REAL*8, ALLOCATABLE :: a_tilde_ij(:,:)
130  REAL*8, ALLOCATABLE :: a_dirk_ij(:,:)
131 
133  REAL*8, ALLOCATABLE :: omega_tilde(:)
134 
136  REAL*8, ALLOCATABLE :: omega(:)
137 
139  REAL*8, ALLOCATABLE :: a_tilde(:)
140 
142  REAL*8, ALLOCATABLE :: a_dirk(:)
143 
145  REAL*8 :: a_diag
146 
148  REAL*8, ALLOCATABLE :: q_rk(:,:,:,:)
149 
151  REAL*8, ALLOCATABLE :: qp_rk(:,:,:,:)
152 
153 
155  REAL*8, ALLOCATABLE :: divflux(:,:,:,:)
156 
158  REAL*8, ALLOCATABLE :: nh(:,:,:,:)
159 
161  REAL*8, ALLOCATABLE :: si_nh(:,:,:,:)
162 
164  REAL*8, ALLOCATABLE :: expl_terms(:,:,:,:)
165 
167  REAL*8, ALLOCATABLE :: divfluxj(:,:)
168 
170  REAL*8, ALLOCATABLE :: nhj(:,:)
171 
173  REAL*8, ALLOCATABLE :: expl_terms_j(:,:)
174 
176  REAL*8, ALLOCATABLE :: si_nhj(:,:)
177 
179  LOGICAL :: normalize_q
180 
182  LOGICAL :: normalize_f
183 
185  LOGICAL :: opt_search_nl
186 
188  REAL*8, ALLOCATABLE :: residual_term(:,:,:)
189 
190  INTEGER, ALLOCATABLE :: j_cent(:)
191  INTEGER, ALLOCATABLE :: k_cent(:)
192 
193  INTEGER, ALLOCATABLE :: j_stag_x(:)
194  INTEGER, ALLOCATABLE :: k_stag_x(:)
195 
196  INTEGER, ALLOCATABLE :: j_stag_y(:)
197  INTEGER, ALLOCATABLE :: k_stag_y(:)
198 
199  REAL*8 :: t_imex1,t_imex2
200 
201 CONTAINS
202 
203  !******************************************************************************
205  !
208  !
212  !
213  !******************************************************************************
214 
215  SUBROUTINE allocate_solver_variables
216 
217  IMPLICIT NONE
218 
219  REAL*8 :: gamma , delta
220 
221  INTEGER :: i,j
222 
223  ALLOCATE( q( n_vars , comp_cells_x , comp_cells_y ) , q0( n_vars , &
224  comp_cells_x , comp_cells_y ) )
225 
226  ALLOCATE( qp( n_vars+2 , comp_cells_x , comp_cells_y ) )
227 
228  q(1:n_vars,1:comp_cells_x,1:comp_cells_y) = 0.d0
229  qp(1:n_vars+2,1:comp_cells_x,1:comp_cells_y) = 0.d0
230 
231  ALLOCATE( q1max( comp_cells_x , comp_cells_y ) )
232 
233  ALLOCATE( q_fv( n_vars , comp_cells_x , comp_cells_y ) )
234 
235  ALLOCATE( q_interfacel( n_vars , comp_interfaces_x, comp_cells_y ) )
236  ALLOCATE( q_interfacer( n_vars , comp_interfaces_x, comp_cells_y ) )
237  ALLOCATE( a_interface_xneg( n_eqns , comp_interfaces_x, comp_cells_y ) )
238  ALLOCATE( a_interface_xpos( n_eqns , comp_interfaces_x, comp_cells_y ) )
239  ALLOCATE( h_interface_x( n_eqns , comp_interfaces_x, comp_cells_y ) )
240 
241 
242  ALLOCATE( q_interfaceb( n_vars , comp_cells_x, comp_interfaces_y ) )
243  ALLOCATE( q_interfacet( n_vars , comp_cells_x, comp_interfaces_y ) )
244  ALLOCATE( a_interface_yneg( n_eqns , comp_cells_x, comp_interfaces_y ) )
245  ALLOCATE( a_interface_ypos( n_eqns , comp_cells_x, comp_interfaces_y ) )
246 
247  ALLOCATE( qp_interfacel( n_vars+2 , comp_interfaces_x, comp_cells_y ) )
248  ALLOCATE( qp_interfacer( n_vars+2 , comp_interfaces_x, comp_cells_y ) )
249  ALLOCATE( qp_interfaceb( n_vars+2 , comp_cells_x, comp_interfaces_y ) )
250  ALLOCATE( qp_interfacet( n_vars+2 , comp_cells_x, comp_interfaces_y ) )
251 
252 
253  ALLOCATE( q_cellnw( n_vars , comp_cells_x , comp_cells_y ) )
254  ALLOCATE( q_cellne( n_vars , comp_cells_x , comp_cells_y ) )
255  ALLOCATE( q_cellsw( n_vars , comp_cells_x , comp_cells_y ) )
256  ALLOCATE( q_cellse( n_vars , comp_cells_x , comp_cells_y ) )
257 
258  ALLOCATE( qp_cellnw( n_vars+2 , comp_cells_x , comp_cells_y ) )
259  ALLOCATE( qp_cellne( n_vars+2 , comp_cells_x , comp_cells_y ) )
260  ALLOCATE( qp_cellsw( n_vars+2 , comp_cells_x , comp_cells_y ) )
261  ALLOCATE( qp_cellse( n_vars+2 , comp_cells_x , comp_cells_y ) )
262 
263 
264  ALLOCATE ( a_interface_x_max(n_eqns,comp_interfaces_x,comp_cells_y) )
265  ALLOCATE ( a_interface_y_max(n_eqns,comp_cells_x,comp_interfaces_y) )
266 
267 
268  ALLOCATE( h_interface_y( n_eqns , comp_cells_x, comp_interfaces_y ) )
269 
270  ALLOCATE( solve_mask( comp_cells_x , comp_cells_y ) )
271  ALLOCATE( solve_mask_x( comp_interfaces_x , comp_cells_y ) )
272  ALLOCATE( solve_mask_y( comp_cells_x , comp_interfaces_y ) )
273 
274  ALLOCATE( source_xy( comp_cells_x , comp_cells_y ) )
275 
276 
277  ALLOCATE( a_tilde_ij(n_rk,n_rk) )
278  ALLOCATE( a_dirk_ij(n_rk,n_rk) )
279  ALLOCATE( omega_tilde(n_rk) )
280  ALLOCATE( omega(n_rk) )
281 
282 
283  ! Allocate the logical arrays defining the implicit parts of the system
284  ALLOCATE( mask22(n_eqns,n_eqns) )
285  ALLOCATE( mask21(n_eqns,n_eqns) )
286  ALLOCATE( mask11(n_eqns,n_eqns) )
287  ALLOCATE( mask12(n_eqns,n_eqns) )
288 
289  ! Initialize the logical arrays with all false (everything is implicit)
290  mask11(1:n_eqns,1:n_eqns) = .false.
291  mask12(1:n_eqns,1:n_eqns) = .false.
292  mask22(1:n_eqns,1:n_eqns) = .false.
293  mask21(1:n_eqns,1:n_eqns) = .false.
294 
295  ! Set to .TRUE. the elements not corresponding to equations and variables to
296  ! be solved implicitly
297  DO i = 1,n_eqns
298 
299  DO j = 1,n_eqns
300 
301  IF ( .NOT.implicit_flag(i) .AND. .NOT.implicit_flag(j) ) &
302  mask11(j,i) = .true.
303  IF ( implicit_flag(i) .AND. .NOT.implicit_flag(j) ) &
304  mask12(j,i) = .true.
305  IF ( implicit_flag(i) .AND. implicit_flag(j) ) &
306  mask22(j,i) = .true.
307  IF ( .NOT.implicit_flag(i) .AND. implicit_flag(j) ) &
308  mask21(j,i) = .true.
309 
310  END DO
311 
312  END DO
313 
314  ! Initialize the coefficients for the IMEX Runge-Kutta scheme
315  ! Please note that with respect to the schemes described in Pareschi & Russo
316  ! (2000) we do not have the coefficient vectors c_tilde and c, because the
317  ! explicit and implicit terms do not depend explicitly on time.
318 
319  ! Explicit part coefficients (a_tilde_ij=0 for j>=i)
320  a_tilde_ij = 0.d0
321 
322  ! Weight coefficients of the explicit part in the final assemblage
323  omega_tilde = 0.d0
324 
325  ! Implicit part coefficients (a_dirk_ij=0 for j>i)
326  a_dirk_ij = 0.d0
327 
328  ! Weight coefficients of the explicit part in the final assemblage
329  omega = 0.d0
330 
331  gamma = 1.d0 - 1.d0 / sqrt(2.d0)
332  delta = 1.d0 - 1.d0 / ( 2.d0 * gamma )
333 
334  IF ( n_rk .EQ. 1 ) THEN
335 
336  a_tilde_ij(1,1) = 1.d0
337 
338  omega_tilde(1) = 1.d0
339 
340  a_dirk_ij(1,1) = 0.d0
341 
342  omega(1) = 0.d0
343 
344  ELSEIF ( n_rk .EQ. 2 ) THEN
345 
346  a_tilde_ij(2,1) = 1.0d0
347 
348  omega_tilde(1) = 1.0d0
349  omega_tilde(2) = 0.0d0
350 
351  a_dirk_ij(2,2) = 1.0d0
352 
353  omega(1) = 0.d0
354  omega(2) = 1.d0
355 
356 
357  ELSEIF ( n_rk .EQ. 3 ) THEN
358 
359  ! Tableau for the IMEX-SSP(3,3,2) Stiffly Accurate Scheme
360  ! from Pareschi & Russo (2005), Table IV
361 
362  a_tilde_ij(2,1) = 0.5d0
363  a_tilde_ij(3,1) = 0.5d0
364  a_tilde_ij(3,2) = 0.5d0
365 
366  omega_tilde(1) = 1.0d0 / 3.0d0
367  omega_tilde(2) = 1.0d0 / 3.0d0
368  omega_tilde(3) = 1.0d0 / 3.0d0
369 
370  a_dirk_ij(1,1) = 0.25d0
371  a_dirk_ij(2,2) = 0.25d0
372  a_dirk_ij(3,1) = 1.0d0 / 3.0d0
373  a_dirk_ij(3,2) = 1.0d0 / 3.0d0
374  a_dirk_ij(3,3) = 1.0d0 / 3.0d0
375 
376  omega(1) = 1.0d0 / 3.0d0
377  omega(2) = 1.0d0 / 3.0d0
378  omega(3) = 1.0d0 / 3.0d0
379 
380  ELSEIF ( n_rk .EQ. 4 ) THEN
381 
382  ! LRR(3,2,2) from Table 3 in Pareschi & Russo (2000)
383 
384  a_tilde_ij(2,1) = 0.5d0
385  a_tilde_ij(3,1) = 1.d0 / 3.d0
386  a_tilde_ij(4,2) = 1.0d0
387 
388  omega_tilde(1) = 0.d0
389  omega_tilde(2) = 1.0d0
390  omega_tilde(3) = 0.0d0
391  omega_tilde(4) = 0.d0
392 
393  a_dirk_ij(2,2) = 0.5d0
394  a_dirk_ij(3,3) = 1.0d0 / 3.0d0
395  a_dirk_ij(4,3) = 0.75d0
396  a_dirk_ij(4,4) = 0.25d0
397 
398  omega(1) = 0.d0
399  omega(2) = 0.d0
400  omega(3) = 0.75d0
401  omega(4) = 0.25d0
402 
403  END IF
404 
405  ALLOCATE( a_tilde(n_rk) )
406  ALLOCATE( a_dirk(n_rk) )
407 
408  ALLOCATE( q_rk( n_vars , comp_cells_x , comp_cells_y , n_rk ) )
409  ALLOCATE( qp_rk( n_vars+2 , comp_cells_x , comp_cells_y , n_rk ) )
410  ALLOCATE( divflux( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )
411  ALLOCATE( nh( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )
412  ALLOCATE( si_nh( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )
413 
414  ALLOCATE( expl_terms( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )
415 
416  ALLOCATE( divfluxj(n_eqns,n_rk) )
417  ALLOCATE( nhj(n_eqns,n_rk) )
418  ALLOCATE( si_nhj(n_eqns,n_rk) )
419 
420  ALLOCATE( expl_terms_j(n_eqns,n_rk) )
421 
422  ALLOCATE( residual_term( n_vars , comp_cells_x , comp_cells_y ) )
423 
424  comp_cells_xy = comp_cells_x * comp_cells_y
425 
426  ALLOCATE( j_cent( comp_cells_xy ) )
427  ALLOCATE( k_cent( comp_cells_xy ) )
428 
429  ALLOCATE( j_stag_x( comp_interfaces_x * comp_cells_y ) )
430  ALLOCATE( k_stag_x( comp_interfaces_x * comp_cells_y ) )
431 
432  ALLOCATE( j_stag_y( comp_cells_x * comp_interfaces_y ) )
433  ALLOCATE( k_stag_y( comp_cells_x * comp_interfaces_y ) )
434 
435  RETURN
436 
437  END SUBROUTINE allocate_solver_variables
438 
439  !******************************************************************************
441  !
444  !
448  !
449  !******************************************************************************
450 
451  SUBROUTINE deallocate_solver_variables
452 
453  DEALLOCATE( q , q0 )
454 
455  DEALLOCATE( q1max )
456 
457  DEALLOCATE( q_fv )
458 
459  DEALLOCATE( q_interfacel )
460  DEALLOCATE( q_interfacer )
461  DEALLOCATE( q_interfaceb )
462  DEALLOCATE( q_interfacet )
463 
464  DEALLOCATE( q_cellnw )
465  DEALLOCATE( q_cellne )
466  DEALLOCATE( q_cellsw )
467  DEALLOCATE( q_cellse )
468 
469  DEALLOCATE( qp_interfacel )
470  DEALLOCATE( qp_interfacer )
471  DEALLOCATE( qp_interfaceb )
472  DEALLOCATE( qp_interfacet )
473 
474  DEALLOCATE( qp_cellnw )
475  DEALLOCATE( qp_cellne )
476  DEALLOCATE( qp_cellsw )
477  DEALLOCATE( qp_cellse )
478 
479 
480  DEALLOCATE( a_interface_xneg )
481  DEALLOCATE( a_interface_xpos )
482  DEALLOCATE( a_interface_yneg )
483  DEALLOCATE( a_interface_ypos )
484 
485  DEALLOCATE( a_interface_x_max )
486  DEALLOCATE( a_interface_y_max )
487 
488  DEALLOCATE( h_interface_x )
489  DEALLOCATE( h_interface_y )
490 
491  DEALLOCATE( solve_mask )
492  DEALLOCATE( solve_mask_x )
493  DEALLOCATE( solve_mask_y )
494 
495  Deallocate( qp )
496 
497  DEALLOCATE( source_xy )
498 
499  DEALLOCATE( a_tilde_ij )
500  DEALLOCATE( a_dirk_ij )
501  DEALLOCATE( omega_tilde )
502  DEALLOCATE( omega )
503 
504  DEALLOCATE( implicit_flag )
505 
506  DEALLOCATE( a_tilde )
507  DEALLOCATE( a_dirk )
508 
509  DEALLOCATE( q_rk )
510  DEALLOCATE( qp_rk )
511  DEALLOCATE( divflux )
512  DEALLOCATE( nh )
513  DEALLOCATE( si_nh )
514  DEALLOCATE( expl_terms )
515 
516  DEALLOCATE( divfluxj )
517  DEALLOCATE( nhj )
518  DEALLOCATE( si_nhj )
519  DEALLOCATE( expl_terms_j )
520 
521  DEALLOCATE( mask22 , mask21 , mask11 , mask12 )
522 
523  DEALLOCATE( residual_term )
524 
525  DEALLOCATE( j_cent , k_cent )
526  DEALLOCATE ( j_stag_x , k_stag_x )
527  DEALLOCATE ( j_stag_y , k_stag_y )
528 
529  RETURN
530 
531  END SUBROUTINE deallocate_solver_variables
532 
533 
534  !******************************************************************************
536  !
540  !
544  !
545  !******************************************************************************
546 
547  SUBROUTINE check_solve
548 
549  IMPLICIT NONE
550 
551  INTEGER :: i,j,k
552 
553  solve_mask(1:comp_cells_x,1:comp_cells_y) = .false.
554 
555  WHERE ( q(1,:,:) .GT. 0.d0 ) solve_mask = .true.
556 
557  ! the equations are solved here to prescribe boundary conditions
558  solve_mask(1,1:comp_cells_y) = .true.
559  solve_mask(comp_cells_x,1:comp_cells_y) = .true.
560  solve_mask(1:comp_cells_x,1) = .true.
561  solve_mask(1:comp_cells_x,comp_cells_y) = .true.
562 
563  IF ( radial_source_flag ) THEN
564 
565  DO k = 1,comp_cells_y
566 
567  DO j = 1,comp_cells_x
568 
569  IF ( source_cell(j,k) .EQ. 2 ) THEN
570 
571  solve_mask(j,k) = .true.
572 
573  END IF
574 
575  END DO
576 
577  END DO
578 
579  END IF
580 
581  DO i = 1,n_rk
582 
583  ! solution domain is extended to neighbours of positive-mass cells
584  solve_mask(2:comp_cells_x-1,2:comp_cells_y-1) = &
585  solve_mask(2:comp_cells_x-1,2:comp_cells_y-1) .OR. &
586  solve_mask(1:comp_cells_x-2,2:comp_cells_y-1) .OR. &
587  solve_mask(3:comp_cells_x,2:comp_cells_y-1) .OR. &
588  solve_mask(2:comp_cells_x-1,1:comp_cells_y-2) .OR. &
589  solve_mask(2:comp_cells_x-1,3:comp_cells_y)
590 
591  END DO
592 
593  IF ( radial_source_flag ) THEN
594 
595  DO k = 1,comp_cells_y
596 
597  DO j = 1,comp_cells_x
598 
599  IF ( source_cell(j,k) .EQ. 1 ) solve_mask(j,k) = .false.
600 
601  END DO
602 
603  END DO
604 
605  END IF
606 
607  solve_mask_x(1:comp_interfaces_x,1:comp_cells_y) = .false.
608  solve_mask_y(1:comp_cells_x,1:comp_interfaces_y) = .false.
609 
610  !----- check for cells where computation is needed
611  i = 0
612 
613  DO k = 1,comp_cells_y
614 
615  DO j = 1,comp_cells_x
616 
617  IF ( solve_mask(j,k) ) THEN
618 
619  i = i+1
620  j_cent(i) = j
621  k_cent(i) = k
622 
623  solve_mask_x(j,k) = .true.
624  solve_mask_x(j+1,k) = .true.
625  solve_mask_y(j,k) = .true.
626  solve_mask_y(j,k+1) = .true.
627 
628  END IF
629 
630  END DO
631 
632  END DO
633 
634  solve_cells = i
635 
636  !----- check for y-interfaces where computation is needed
637  i = 0
638 
639  DO k = 1,comp_cells_y
640 
641  DO j = 1,comp_interfaces_x
642 
643  IF ( solve_mask_x(j,k) ) THEN
644 
645  i = i+1
646  j_stag_x(i) = j
647  k_stag_x(i) = k
648 
649  END IF
650 
651  END DO
652 
653  END DO
654 
655  solve_interfaces_x = i
656 
657  !----- check for y-interfaces where computation is needed
658 
659  i = 0
660 
661  DO k = 1,comp_interfaces_y
662 
663  DO j = 1,comp_cells_x
664 
665  IF ( solve_mask_y(j,k) ) THEN
666 
667  i = i+1
668  j_stag_y(i) = j
669  k_stag_y(i) = k
670 
671  END IF
672 
673  END DO
674 
675  END DO
676 
677  solve_interfaces_y = i
678 
679  RETURN
680 
681  END SUBROUTINE check_solve
682 
683  !*****************************************************************************
685  !
689  !
693  !
694  !*****************************************************************************
695 
696  SUBROUTINE timestep
697 
698  ! External variables
699  USE geometry_2d, ONLY : dx,dy
700  USE parameters_2d, ONLY : max_dt , cfl
701 
702  USE constitutive_2d, ONLY : qc_to_qp
703 
704  IMPLICIT none
705 
706  REAL*8 :: dt_cfl
707 
708  REAL*8 :: dt_interface_x, dt_interface_y
709 
710  INTEGER :: i,j,k,l
711 
712  REAL*8 :: max_a
713 
714  REAL*8 :: max_vel
715  INTEGER :: max_j,max_k
716  INTEGER :: max_dir
717 
718  max_vel = 0.d0
719 
720  dt = max_dt
721 
722  IF ( cfl .NE. -1.d0 ) THEN
723 
724  DO l = 1,solve_cells
725 
726  j = j_cent(l)
727  k = k_cent(l)
728 
729  CALL qc_to_qp( q(1:n_vars,j,k) , qp(1:n_vars+2,j,k) )
730 
731  END DO
732 
733 
734  ! Compute the physical and conservative variables at the interfaces
735  CALL reconstruction( q , qp )
736 
737  ! Compute the max/min eigenvalues at the interfaces
738  CALL eval_speeds
739 
740  DO i=1,n_vars
741 
742  a_interface_x_max(i,:,:) = &
743  max( a_interface_xpos(i,:,:) , -a_interface_xneg(i,:,:) )
744 
745  a_interface_y_max(i,:,:) = &
746  max( a_interface_ypos(i,:,:) , -a_interface_yneg(i,:,:) )
747 
748  END DO
749 
750  DO l = 1,solve_cells
751 
752  j = j_cent(l)
753  k = k_cent(l)
754 
755 
756  max_a = max( maxval(a_interface_x_max(:,j,k)) , &
757  maxval(a_interface_x_max(:,j+1,k)) )
758 
759  IF ( max_a .GT. max_vel) THEN
760 
761  max_vel = max_a
762  max_j = j
763  max_k = k
764  max_dir = 1
765 
766  END IF
767 
768  IF ( max_a .GT. 0.d0 ) THEN
769 
770  dt_interface_x = cfl * dx / max_a
771 
772  ELSE
773 
774  dt_interface_x = dt
775 
776  END IF
777 
778  max_a = max( maxval(a_interface_y_max(:,j,k)) , &
779  maxval(a_interface_y_max(:,j,k+1)) )
780 
781  IF ( max_a .GT. max_vel) THEN
782 
783  max_vel = max_a
784  max_j = j
785  max_k = k
786  max_dir = 2
787 
788  END IF
789 
790  IF ( max_a .GT. 0.d0 ) THEN
791 
792  dt_interface_y = cfl * dy / max_a
793 
794  ELSE
795 
796  dt_interface_y = dt
797 
798  END IF
799 
800 
801  dt_cfl = min( dt_interface_x , dt_interface_y )
802 
803  dt = min(dt,dt_cfl)
804 
805  END DO
806 
807  END IF
808 
809  RETURN
810 
811  END SUBROUTINE timestep
812 
813  !******************************************************************************
815  !
820  !
824  !
825  !******************************************************************************
826 
827  SUBROUTINE imex_rk_solver
828 
829  USE constitutive_2d, ONLY : eval_nonhyperbolic_terms
830 
831  USE constitutive_2d, ONLY : eval_nh_semi_impl_terms
832 
833  USE constitutive_2d, ONLY : qc_to_qp
834 
835  IMPLICIT NONE
836 
837  REAL*8 :: q_si(n_vars)
838  REAL*8 :: q_guess(n_vars)
839  INTEGER :: j,k,l
840 
841  REAL*8 :: h_new
842 
843  ! Initialization of the solution guess
844  q0( 1:n_vars , 1:comp_cells_x , 1:comp_cells_y ) = &
845  q( 1:n_vars , 1:comp_cells_x , 1:comp_cells_y )
846 
847  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solver, imex_RK_solver: beginning'
848 
849  ! Initialization of the variables for the Runge-Kutta scheme
850  q_rk(1:n_vars,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0
851  qp_rk(1:n_vars+2,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0
852 
853  divflux(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0
854 
855  nh(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0
856 
857  si_nh(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0
858 
859  expl_terms(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0
860 
861  runge_kutta:DO i_rk = 1,n_rk
862 
863  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solver, imex_RK_solver: i_RK',i_rk
864 
865  ! define the explicits coefficients for the i-th step of the Runge-Kutta
866  a_tilde = 0.d0
867  a_dirk = 0.d0
868 
869  ! in the first step of the RK scheme all the coefficients remain to 0
870  a_tilde(1:i_rk-1) = a_tilde_ij(i_rk,1:i_rk-1)
871  a_dirk(1:i_rk-1) = a_dirk_ij(i_rk,1:i_rk-1)
872 
873  ! define the implicit coefficient for the i-th step of the Runge-Kutta
874  a_diag = a_dirk_ij(i_rk,i_rk)
875 
876  CALL cpu_time(t_imex1)
877 
878  DO l = 1,solve_cells
879 
880  j = j_cent(l)
881  k = k_cent(l)
882 
883  IF ( verbose_level .GE. 2 ) THEN
884 
885  WRITE(*,*) 'solver, imex_RK_solver: j',j,k
886 
887  END IF
888 
889  ! initialize the RK step
890  IF ( i_rk .EQ. 1 ) THEN
891 
892  ! solution from the previous time step
893  q_guess(1:n_vars) = q0( 1:n_vars , j , k)
894 
895  ELSE
896 
897  ! solution from the previous RK step
898  !q_guess(1:n_vars) = q_rk( 1:n_vars , j , k , MAX(1,i_RK-1) )
899 
900  END IF
901 
902  ! RK explicit hyperbolic terms for the volume (i,j)
903  divfluxj(1:n_eqns,1:n_rk) = divflux( 1:n_eqns , j , k , 1:n_rk )
904 
905  ! RK implicit terms for the volume (i,j)
906  nhj(1:n_eqns,1:n_rk) = nh( 1:n_eqns , j , k , 1:n_rk )
907 
908  ! RK semi-implicit terms for the volume (i,j)
909  si_nhj(1:n_eqns,1:n_rk) = si_nh( 1:n_eqns , j , k , 1:n_rk )
910 
911  ! RK additional explicit terms for the volume (i,j)
912  expl_terms_j(1:n_eqns,1:n_rk) = expl_terms( 1:n_eqns,j,k,1:n_rk )
913 
914  ! New solution at the i_RK step without the implicit and
915  ! semi-implicit term
916  q_fv( 1:n_vars , j , k ) = q0( 1:n_vars , j , k ) &
917  - dt * (matmul( divfluxj(1:n_eqns,1:i_rk) &
918  + expl_terms_j(1:n_eqns,1:i_rk) , a_tilde(1:i_rk) ) &
919  - matmul( nhj(1:n_eqns,1:i_rk) + si_nhj(1:n_eqns,1:i_rk) , &
920  a_dirk(1:i_rk) ) )
921 
922  IF ( verbose_level .GE. 2 ) THEN
923 
924  WRITE(*,*) 'q_guess',q_guess
925  CALL qc_to_qp( q_guess , qp(1:n_vars+2,j,k) )
926  WRITE(*,*) 'q_guess: qp',qp(1:n_vars+2,j,k)
927 
928  END IF
929 
930  adiag_pos:IF ( a_diag .NE. 0.d0 ) THEN
931 
932  pos_thick:IF ( q_fv(1,j,k) .GT. 0.d0 ) THEN
933 
934  ! Eval the semi-implicit terms
935  ! (terms which non depend on velocity magnitude)
936  CALL eval_nh_semi_impl_terms( grav_surf(j,k) , &
937  q_fv( 1:n_vars , j , k ) , si_nh(1:n_eqns,j,k,i_rk) )
938 
939  si_nhj(1:n_eqns,i_rk) = si_nh( 1:n_eqns,j,k,i_rk )
940 
941  ! Assemble the initial guess for the implicit solver
942  q_si(1:n_vars) = q_fv(1:n_vars,j,k ) + dt * a_diag * &
943  si_nh(1:n_eqns,j,k,i_rk)
944 
945  IF ( ( q_fv(2,j,k)**2 + q_fv(3,j,k)**2 ) .EQ. 0.d0 ) THEN
946 
947  !Case 1: if the velocity was null, then it must stay null
948  q_si(2:3) = 0.d0
949 
950  ELSEIF ( ( q_si(2)*q_fv(2,j,k) .LT. 0.d0 ) .OR. &
951  ( q_si(3)*q_fv(3,j,k) .LT. 0.d0 ) ) THEN
952 
953  ! If the semi-impl. friction term changed the sign of the
954  ! velocity then set it to zero
955  q_si(2:3) = 0.d0
956 
957  ELSE
958 
959  ! Align the velocity vector with previous one
960  q_si(2:3) = dsqrt( q_si(2)**2 + q_si(3)**2 ) * &
961  q_fv(2:3,j,k) / dsqrt( q_fv(2,j,k)**2 &
962  + q_fv(3,j,k)**2 )
963 
964  END IF
965 
966  ! Update the semi-implicit term accordingly with the
967  ! corrections above
968  si_nh(1:n_eqns,j,k,i_rk) = ( q_si(1:n_vars) - &
969  q_fv(1:n_vars,j,k ) ) / ( dt*a_diag )
970 
971  si_nhj(1:n_eqns,i_rk) = si_nh( 1:n_eqns,j,k,i_rk )
972 
973  ! Initialize the guess for the NR solver
974  q_guess(1:n_vars) = q_si(1:n_vars)
975 
976  ! Solve the implicit system to find the solution at the
977  ! i_RK step of the IMEX RK procedure
978  CALL solve_rk_step( q_guess(1:n_vars) , q0(1:n_vars,j,k ) , &
979  a_tilde , a_dirk , a_diag )
980 
981  IF ( comp_cells_y .EQ. 1 ) THEN
982 
983  q_guess(3) = 0.d0
984 
985  END IF
986 
987  IF ( comp_cells_x .EQ. 1 ) THEN
988 
989  q_guess(2) = 0.d0
990 
991  END IF
992 
993  ! Eval and store the implicit term at the i_RK step
994  CALL eval_nonhyperbolic_terms( r_qj =q_guess , &
995  r_nh_term_impl = nh(1:n_eqns,j,k,i_rk) )
996 
997  IF ( q_si(2)**2 + q_si(3)**2 .EQ. 0.d0 ) THEN
998 
999  q_guess(2:3) = 0.d0
1000 
1001  ELSEIF ( ( q_guess(2)*q_si(2) .LE. 0.d0 ) .AND. &
1002  ( q_guess(3)*q_si(3) .LE. 0.d0 ) ) THEN
1003 
1004  ! If the impl. friction term changed the sign of the
1005  ! velocity then set it to zero
1006  q_guess(2:3) = 0.d0
1007 
1008  ELSE
1009 
1010  ! Align the velocity vector with previous one
1011  q_guess(2:3) = dsqrt( q_guess(2)**2 + q_guess(3)**2 ) * &
1012  q_si(2:3) / dsqrt( q_si(2)**2 + q_si(3)**2 )
1013 
1014  END IF
1015 
1016  ELSE
1017 
1018  ! If h=0 nothing has to be changed
1019  q_guess(1:n_vars) = q_fv( 1:n_vars , j , k )
1020  q_si(1:n_vars) = q_fv( 1:n_vars , j , k )
1021  si_nh(1:n_eqns,j,k,i_rk) = 0.d0
1022  nh(1:n_eqns,j,k,i_rk) = 0.d0
1023 
1024  END IF pos_thick
1025 
1026  END IF adiag_pos
1027 
1028  ! Check the sign of the flow thickness
1029  h_new = q_guess(1)
1030 
1031  IF ( a_diag .NE. 0.d0 ) THEN
1032 
1033  ! Update the implicit term with correction on the new velocity
1034  nh(1:n_vars,j,k,i_rk) = ( q_guess(1:n_vars) - q_si(1:n_vars)) &
1035  / ( dt*a_diag )
1036 
1037  END IF
1038 
1039  ! Store the solution at the end of the i_RK step
1040  q_rk( 1:n_vars , j , k , i_rk ) = q_guess
1041 
1042  IF ( verbose_level .GE. 2 ) THEN
1043 
1044  WRITE(*,*) 'imex_RK_solver: qc',q_guess
1045  CALL qc_to_qp( q_guess , qp(1:n_vars+2,j,k) )
1046  WRITE(*,*) 'imex_RK_solver: qp',qp(1:n_vars+2,j,k)
1047  READ(*,*)
1048 
1049  END IF
1050 
1051  END DO
1052 
1053  CALL cpu_time(t_imex2)
1054  ! WRITE(*,*) 'Time taken by implicit',t_imex2-t_imex1,'seconds'
1055 
1056  IF ( omega_tilde(i_rk) .GT. 0.d0 ) THEN
1057 
1058  ! Compute the physical variables at the cell centers
1059  DO l = 1,solve_cells
1060 
1061  j = j_cent(l)
1062  k = k_cent(l)
1063 
1064  CALL qc_to_qp( q_rk(1:n_vars,j,k,i_rk) , &
1065  qp_rk(1:n_vars+2,j,k,i_rk) )
1066 
1067  END DO
1068 
1069 
1070  ! Eval and store the explicit hyperbolic (fluxes) terms
1071  CALL eval_hyperbolic_terms( &
1072  q_rk(1:n_vars,1:comp_cells_x,1:comp_cells_y,i_rk) , &
1073  qp_rk(1:n_vars+2,1:comp_cells_x,1:comp_cells_y,i_rk) , &
1074  divflux(1:n_eqns,1:comp_cells_x,1:comp_cells_y,i_rk) )
1075 
1076  CALL cpu_time(t_imex1)
1077  ! WRITE(*,*) 'Time taken by explicit',t_imex1-t_imex2,'seconds'
1078 
1079  ! Eval and store the other explicit terms (e.g. gravity or viscous
1080  ! forces)
1081  CALL eval_explicit_terms( &
1082  q_rk(1:n_vars,1:comp_cells_x,1:comp_cells_y,i_rk) , &
1083  qp_rk(1:n_vars+2,1:comp_cells_x,1:comp_cells_y,i_rk) , &
1084  expl_terms(1:n_eqns,1:comp_cells_x,1:comp_cells_y,i_rk) )
1085 
1086  CALL cpu_time(t_imex1)
1087  ! WRITE(*,*) 'Time taken by explicit',t_imex1-t_imex2,'seconds'
1088 
1089  END IF
1090 
1091  END DO runge_kutta
1092 
1093  DO l = 1,solve_cells
1094 
1095  j = j_cent(l)
1096  k = k_cent(l)
1097 
1098  residual_term(1:n_vars,j,k) = matmul( divflux(1:n_eqns,j,k,1:n_rk) &
1099  + expl_terms(1:n_eqns,j,k,1:n_rk) , omega_tilde ) - &
1100  matmul( nh(1:n_eqns,j,k,1:n_rk) + si_nh(1:n_eqns,j,k,1:n_rk) , &
1101  omega )
1102 
1103  END DO
1104 
1105  assemble_sol:DO l = 1,solve_cells
1106 
1107  j = j_cent(l)
1108  k = k_cent(l)
1109 
1110  IF ( verbose_level .GE. 1 ) THEN
1111 
1112  WRITE(*,*) 'cell jk =',j,k
1113  WRITE(*,*) 'before imex_RK_solver: qc',q0(1:n_vars,j,k)
1114  CALL qc_to_qp(q0(1:n_vars,j,k) , qp(1:n_vars+2,j,k))
1115  WRITE(*,*) 'before imex_RK_solver: qp',qp(1:n_vars+2,j,k)
1116 
1117  END IF
1118 
1119  IF ( ( sum(abs( omega_tilde(:)-a_tilde_ij(n_rk,:))) .EQ. 0.d0 ) &
1120  .AND. ( sum(abs(omega(:)-a_dirk_ij(n_rk,:))) .EQ. 0.d0 ) ) THEN
1121 
1122  ! The assembling coeffs are equal to the last step of the RK scheme
1123  q(1:n_vars,j,k) = q_rk(1:n_vars,j,k,n_rk)
1124 
1125  ELSE
1126 
1127  ! The assembling coeffs are different
1128  q(1:n_vars,j,k) = q0(1:n_vars,j,k) - dt*residual_term(1:n_vars,j,k)
1129 
1130  END IF
1131 
1132  negative_thickness_check:IF ( q(1,j,k) .LT. 0.d0 ) THEN
1133 
1134  IF ( q(1,j,k) .GT. -1.d-7 ) THEN
1135 
1136  q(1,j,k) = 0.d0
1137  q(2:n_vars,j,k) = 0.d0
1138 
1139  ELSE
1140 
1141  WRITE(*,*) 'j,k,n_RK',j,k,n_rk
1142  WRITE(*,*) 'dt',dt
1143  WRITE(*,*) 'source_cell',source_cell(j,k)
1144  WRITE(*,*) 'source_cell left',source_cell(j-1,k)
1145  WRITE(*,*) 'source_cell right',source_cell(j+1,k)
1146  WRITE(*,*) 'source_cell top',source_cell(j,k+1)
1147  WRITE(*,*) 'source_cell bottom',source_cell(j,k-1)
1148  WRITE(*,*) 'before imex_RK_solver: qc',q0(1:n_vars,j,k)
1149  CALL qc_to_qp(q0(1:n_vars,j,k) , qp(1:n_vars+2,j,k))
1150  WRITE(*,*) 'before imex_RK_solver: qp',qp(1:n_vars+2,j,k)
1151  WRITE(*,*) 'h old',q0(1,j,k)
1152  WRITE(*,*) 'h new',q(1,j,k)
1153  WRITE(*,*) 'B_cent(j,k)',b_cent(j,k)
1154  WRITE(*,*) 'qp_interfaceR(1:n_vars+2,j,k)',qp_interfacer(1:n_vars+2,j,k)
1155  WRITE(*,*) 'qp_interfaceL(1:n_vars+2,j+1,k)',qp_interfacel(1:n_vars+2,j+1,k)
1156  WRITE(*,*) 'qp_interfaceT(1:n_vars+2,j,k)',qp_interfacet(1:n_vars+2,j,k)
1157  WRITE(*,*) 'qp_interfaceB(1:n_vars+2,j,k+1)',qp_interfaceb(1:n_vars+2,j,k)
1158 
1159  WRITE(*,*) 'hS',q_interfacet(1,j,k)
1160  WRITE(*,*) 'hE',q_interfacer(1,j,k)
1161 
1162  READ(*,*)
1163 
1164  END IF
1165 
1166  END IF negative_thickness_check
1167 
1168  IF ( sum(q(5:4+n_solid,j,k)) .GT. q(1,j,k) ) THEN
1169 
1170  IF ( sum(q(5:4+n_solid,j,k))-q(1,j,k) .LT. 1.d-10 ) THEN
1171 
1172  q(5:4+n_solid,j,k) = q(5:4+n_solid,j,k) &
1173  / sum(q(5:4+n_solid,j,k)) * q(1,j,k)
1174 
1175  ELSE
1176 
1177  WRITE(*,*) 'j,k,n_RK',j,k,n_rk
1178  WRITE(*,*) 'dt',dt
1179  WRITE(*,*) ' B_cent(j,k)', b_cent(j,k)
1180  WRITE(*,*) 'before imex_RK_solver: qc',q0(1:n_vars,j,k)
1181  CALL qc_to_qp(q0(1:n_vars,j,k) , qp(1:n_vars+2,j,k))
1182  WRITE(*,*) 'before imex_RK_solver: qp',qp(1:n_vars+2,j,k)
1183  WRITE(*,*) 'after imex_RK_solver: qc',q(1:n_vars,j,k)
1184  CALL qc_to_qp(q(1:n_vars,j,k) , qp(1:n_vars+2,j,k))
1185  WRITE(*,*) 'after imex_RK_solver: qp',qp(1:n_vars+2,j,k)
1186  READ(*,*)
1187 
1188  END IF
1189 
1190  IF ( verbose_level .GE. 1 ) THEN
1191 
1192  WRITE(*,*) 'h new',q(1,j,k)
1193  READ(*,*)
1194 
1195  END IF
1196 
1197  END IF
1198 
1199  END DO assemble_sol
1200 
1201  RETURN
1202 
1203  END SUBROUTINE imex_rk_solver
1204 
1205  !******************************************************************************
1207  !
1214  !
1220  !
1224  !
1225  !******************************************************************************
1226 
1227  SUBROUTINE solve_rk_step( qj, qj_old, a_tilde , a_dirk , a_diag )
1228 
1229  USE parameters_2d, ONLY : max_nl_iter , tol_rel , tol_abs
1230 
1231  USE constitutive_2d, ONLY : qc_to_qp
1232 
1233  IMPLICIT NONE
1234 
1235  REAL*8, INTENT(INOUT) :: qj(n_vars)
1236  REAL*8, INTENT(IN) :: qj_old(n_vars)
1237  REAL*8, INTENT(IN) :: a_tilde(n_rk)
1238  REAL*8, INTENT(IN) :: a_dirk(n_rk)
1239  REAL*8, INTENT(IN) :: a_diag
1240 
1241  REAL*8 :: qj_init(n_vars)
1242 
1243  REAL*8 :: qj_org(n_vars) , qj_rel(n_vars)
1244 
1245  REAL*8 :: left_matrix(n_eqns,n_vars)
1246  REAL*8 :: right_term(n_eqns)
1247 
1248  REAL*8 :: scal_f
1249 
1250  REAL*8 :: coeff_f(n_eqns)
1251 
1252  REAL*8 :: qj_rel_NR_old(n_vars)
1253  REAL*8 :: scal_f_old
1254  REAL*8 :: desc_dir(n_vars)
1255  REAL*8 :: grad_f(n_vars)
1256  ! REAL*8 :: mod_desc_dir
1257 
1258  INTEGER :: pivot(n_vars)
1259 
1260  REAL*8 :: left_matrix_small22(n_nh,n_nh)
1261  REAL*8 :: left_matrix_small21(n_eqns-n_nh,n_nh)
1262  REAL*8 :: left_matrix_small11(n_eqns-n_nh,n_vars-n_nh)
1263  REAL*8 :: left_matrix_small12(n_nh,n_vars-n_nh)
1264 
1265  REAL*8 :: desc_dir_small2(n_nh)
1266  INTEGER :: pivot_small2(n_nh)
1267 
1268  REAL*8 :: desc_dir_small1(n_vars-n_nh)
1269 
1270  INTEGER :: ok
1271 
1272  INTEGER :: i
1273  INTEGER :: nl_iter
1274 
1275  REAL*8, PARAMETER :: STPMX=100.d0
1276  REAL*8 :: stpmax
1277  LOGICAL :: check
1278 
1279  REAL*8, PARAMETER :: TOLF=1.d-10 , tolmin=1.d-6
1280  REAL*8 :: TOLX
1281 
1282  REAL*8 :: qpj(n_vars+2)
1283 
1284  REAL*8 :: desc_dir2(n_vars)
1285 
1286  REAL*8 :: desc_dir_temp(n_vars)
1287 
1288  normalize_q = .true.
1289  normalize_f = .false.
1290  opt_search_nl = .true.
1291 
1292  coeff_f(1:n_eqns) = 1.d0
1293 
1294  grad_f(1:n_eqns) = 0.d0
1295 
1296  qj_init = qj
1297 
1298  ! normalize the functions of the nonlinear system
1299  IF ( normalize_f ) THEN
1300 
1301  qj = qj_old - dt * ( matmul(divfluxj+ expl_terms_j,a_tilde) &
1302  - matmul(nhj,a_dirk) )
1303 
1304  CALL eval_f( qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , &
1305  right_term , scal_f )
1306 
1307  IF ( verbose_level .GE. 3 ) THEN
1308 
1309  WRITE(*,*) 'solve_rk_step: non-normalized right_term'
1310  WRITE(*,*) right_term
1311  WRITE(*,*) 'scal_f',scal_f
1312 
1313  END IF
1314 
1315  DO i=1,n_eqns
1316 
1317  IF ( abs(right_term(i)) .GE. 1.d0 ) coeff_f(i) = 1.d0 / right_term(i)
1318 
1319  END DO
1320 
1321  right_term = coeff_f * right_term
1322 
1323  scal_f = 0.5d0 * dot_product( right_term , right_term )
1324 
1325  IF ( verbose_level .GE. 3 ) THEN
1326  WRITE(*,*) 'solve_rk_step: after normalization',scal_f
1327  END IF
1328 
1329  END IF
1330 
1331  !---- normalize the conservative variables ------
1332 
1333  IF ( normalize_q ) THEN
1334 
1335  qj_org = qj
1336 
1337  qj_org = max( abs(qj_org) , 1.d-3 )
1338 
1339  ELSE
1340 
1341  qj_org(1:n_vars) = 1.d0
1342 
1343  END IF
1344 
1345  qj_rel = qj / qj_org
1346 
1347  ! -----------------------------------------------
1348 
1349  newton_raphson_loop:DO nl_iter=1,max_nl_iter
1350 
1351  tolx = epsilon(qj_rel)
1352 
1353  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solve_rk_step: nl_iter',nl_iter
1354 
1355  CALL eval_f( qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , &
1356  right_term , scal_f )
1357 
1358  IF ( verbose_level .GE. 2 ) THEN
1359 
1360  WRITE(*,*) 'solve_rk_step: right_term',right_term
1361 
1362  END IF
1363 
1364  IF ( verbose_level .GE. 2 ) THEN
1365 
1366  WRITE(*,*) 'before_lnsrch: scal_f',scal_f
1367 
1368  END IF
1369 
1370  ! check the residual of the system
1371 
1372  IF ( maxval( abs( right_term(:) ) ) < tolf ) THEN
1373 
1374  IF ( verbose_level .GE. 3 ) WRITE(*,*) '1: check',check
1375  RETURN
1376 
1377  END IF
1378 
1379  IF ( ( normalize_f ) .AND. ( scal_f < 1.d-6 ) ) THEN
1380 
1381  IF ( verbose_level .GE. 3 ) WRITE(*,*) 'check scal_f',check
1382  RETURN
1383 
1384  END IF
1385 
1386  ! ---- evaluate the descent direction ------------------------------------
1387 
1388  IF ( count( implicit_flag ) .EQ. n_eqns ) THEN
1389 
1390  CALL eval_jacobian( qj_rel , qj_org,coeff_f , left_matrix )
1391 
1392  desc_dir_temp = - right_term
1393 
1394  CALL dgesv(n_eqns,1, left_matrix , n_eqns, pivot, desc_dir_temp , &
1395  n_eqns, ok)
1396 
1397  desc_dir = desc_dir_temp
1398 
1399  ELSE
1400 
1401  CALL eval_jacobian( qj_rel , qj_org,coeff_f , left_matrix )
1402 
1403  left_matrix_small11 = reshape(pack(left_matrix, mask11), &
1404  [n_eqns-n_nh,n_eqns-n_nh])
1405 
1406  left_matrix_small12 = reshape(pack(left_matrix, mask12), &
1407  [n_nh,n_eqns-n_nh])
1408 
1409  left_matrix_small22 = reshape(pack(left_matrix, mask22), &
1410  [n_nh,n_nh])
1411 
1412  left_matrix_small21 = reshape(pack(left_matrix, mask21), &
1413  [n_eqns-n_nh,n_nh])
1414 
1415  desc_dir_small1 = pack( right_term, .NOT.implicit_flag )
1416  desc_dir_small2 = pack( right_term , implicit_flag )
1417 
1418  DO i=1,n_vars-n_nh
1419 
1420  desc_dir_small1(i) = desc_dir_small1(i) / left_matrix_small11(i,i)
1421 
1422  END DO
1423 
1424  desc_dir_small2 = desc_dir_small2 - &
1425  matmul( desc_dir_small1 , left_matrix_small21 )
1426 
1427  CALL dgesv(n_nh,1, left_matrix_small22 , n_nh , pivot_small2 , &
1428  desc_dir_small2 , n_nh, ok)
1429 
1430  desc_dir = unpack( - desc_dir_small2 , implicit_flag , 0.0d0 ) &
1431  + unpack( - desc_dir_small1 , .NOT.implicit_flag , 0.0d0 )
1432 
1433  END IF
1434 
1435  !mod_desc_dir = DSQRT( desc_dir(2)**2 + desc_dir(3)**2 )
1436  !
1437  !IF ( qj(2)**2 + qj(3)**2 .GT. 0.D0 ) THEN
1438  !
1439  ! desc_dir(2) = mod_desc_dir * qj(2) / ( qj(2)**2 + qj(3)**2 )
1440  ! desc_dir(3) = mod_desc_dir * qj(3) / ( qj(2)**2 + qj(3)**2 )
1441  !
1442  !END IF
1443 
1444  IF ( verbose_level .GE. 3 ) WRITE(*,*) 'desc_dir',desc_dir
1445 
1446  qj_rel_nr_old = qj_rel
1447  scal_f_old = scal_f
1448 
1449  IF ( ( opt_search_nl ) .AND. ( nl_iter .GT. 1 ) ) THEN
1450  ! Search for the step lambda giving a suffic. decrease in the solution
1451 
1452  stpmax = stpmx * max( sqrt( dot_product(qj_rel,qj_rel) ) , &
1453  dble( SIZE(qj_rel) ) )
1454 
1455  grad_f = matmul( right_term , left_matrix )
1456 
1457  desc_dir2 = desc_dir
1458 
1459  CALL lnsrch( qj_rel_nr_old , qj_org , qj_old , scal_f_old , grad_f , &
1460  desc_dir , coeff_f , qj_rel , scal_f , right_term , stpmax , &
1461  check )
1462 
1463  ELSE
1464 
1465  qj_rel = qj_rel_nr_old + desc_dir
1466 
1467  qj = qj_rel * qj_org
1468 
1469  CALL eval_f( qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , &
1470  right_term , scal_f )
1471 
1472  END IF
1473 
1474  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'after_lnsrch: scal_f',scal_f
1475 
1476  qj = qj_rel * qj_org
1477 
1478  IF ( verbose_level .GE. 3 ) THEN
1479 
1480  WRITE(*,*) 'qj',qj
1481  CALL qc_to_qp( qj , qpj)
1482  WRITE(*,*) 'qp',qpj
1483 
1484  END IF
1485 
1486 
1487  IF ( maxval( abs( right_term(:) ) ) < tolf ) THEN
1488 
1489  IF ( verbose_level .GE. 3 ) WRITE(*,*) '1: check',check
1490  check= .false.
1491  RETURN
1492 
1493  END IF
1494 
1495  IF (check) THEN
1496 
1497  check = ( maxval( abs(grad_f(:)) * max( abs( qj_rel(:) ),1.d0 ) / &
1498  max( scal_f , 0.5d0 * SIZE(qj_rel) ) ) < tolmin )
1499 
1500  IF ( verbose_level .GE. 3 ) WRITE(*,*) '2: check',check
1501  ! RETURN
1502 
1503  END IF
1504 
1505  IF ( maxval( abs( qj_rel(:) - qj_rel_nr_old(:) ) / max( abs( qj_rel(:)) ,&
1506  1.d0 ) ) < tolx ) THEN
1507 
1508  IF ( verbose_level .GE. 3 ) WRITE(*,*) 'check',check
1509  RETURN
1510 
1511  END IF
1512 
1513  END DO newton_raphson_loop
1514 
1515  RETURN
1516 
1517  END SUBROUTINE solve_rk_step
1518 
1519  !******************************************************************************
1521  !
1539  !******************************************************************************
1540 
1541  SUBROUTINE lnsrch( qj_rel_NR_old , qj_org , qj_old , scal_f_old , grad_f , &
1542  desc_dir , coeff_f , qj_rel , scal_f , right_term , stpmax , check )
1543 
1544  IMPLICIT NONE
1545 
1547  REAL*8, DIMENSION(:), INTENT(IN) :: qj_rel_NR_old
1548 
1550  REAL*8, DIMENSION(:), INTENT(IN) :: qj_org
1551 
1553  REAL*8, DIMENSION(:), INTENT(IN) :: qj_old
1554 
1556  REAL*8, DIMENSION(:), INTENT(IN) :: grad_f
1557 
1559  REAL*8, INTENT(IN) :: scal_f_old
1560 
1562  REAL*8, DIMENSION(:), INTENT(INOUT) :: desc_dir
1563 
1564  REAL*8, INTENT(IN) :: stpmax
1565 
1567  REAL*8, DIMENSION(:), INTENT(IN) :: coeff_f
1568 
1570  REAL*8, DIMENSION(:), INTENT(OUT) :: qj_rel
1571 
1573  REAL*8, INTENT(OUT) :: scal_f
1574 
1576  REAL*8, INTENT(OUT) :: right_term(n_eqns)
1577 
1579  LOGICAL, INTENT(OUT) :: check
1580 
1581  REAL*8, PARAMETER :: TOLX=epsilon(qj_rel)
1582 
1583  INTEGER, DIMENSION(1) :: ndum
1584  REAL*8 :: ALF , a,alam,alam2,alamin,b,disc
1585  REAL*8 :: scal_f2
1586  REAL*8 :: desc_dir_abs
1587  REAL*8 :: rhs1 , rhs2 , slope, tmplam
1588 
1589  REAL*8 :: scal_f_min , alam_min
1590 
1591  REAL*8 :: qj(n_vars)
1592 
1593  alf = 1.0d-4
1594 
1595  IF ( size(grad_f) == size(desc_dir) .AND. size(grad_f) == size(qj_rel) &
1596  .AND. size(qj_rel) == size(qj_rel_nr_old) ) THEN
1597 
1598  ndum = size(grad_f)
1599 
1600  ELSE
1601 
1602  WRITE(*,*) 'nrerror: an assert_eq failed with this tag:', 'lnsrch'
1603  stop 'program terminated by assert_eq4'
1604 
1605  END IF
1606 
1607  check = .false.
1608 
1609  desc_dir_abs = sqrt( dot_product(desc_dir,desc_dir) )
1610 
1611  IF ( desc_dir_abs > stpmax ) desc_dir(:) = desc_dir(:) * stpmax/desc_dir_abs
1612 
1613  slope = dot_product(grad_f,desc_dir)
1614 
1615  alamin = tolx / maxval( abs( desc_dir(:))/max( abs(qj_rel_nr_old(:)),1.d0 ) )
1616 
1617  IF ( alamin .EQ. 0.d0) THEN
1618 
1619  qj_rel(:) = qj_rel_nr_old(:)
1620 
1621  RETURN
1622 
1623  END IF
1624 
1625  alam = 1.0d0
1626 
1627  scal_f_min = scal_f_old
1628 
1629  optimal_step_search: DO
1630 
1631  IF ( verbose_level .GE. 4 ) THEN
1632 
1633  WRITE(*,*) 'alam',alam
1634 
1635  END IF
1636 
1637  qj_rel = qj_rel_nr_old + alam * desc_dir
1638 
1639  qj = qj_rel * qj_org
1640 
1641  CALL eval_f( qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , &
1642  right_term , scal_f )
1643 
1644  IF ( verbose_level .GE. 4 ) THEN
1645 
1646  WRITE(*,*) 'lnsrch: effe_old,effe',scal_f_old,scal_f
1647  READ(*,*)
1648 
1649  END IF
1650 
1651  IF ( scal_f .LT. scal_f_min ) THEN
1652 
1653  scal_f_min = scal_f
1654  alam_min = alam
1655 
1656  END IF
1657 
1658  IF ( scal_f .LE. 0.9 * scal_f_old ) THEN
1659  ! sufficient function decrease
1660 
1661  IF ( verbose_level .GE. 4 ) THEN
1662 
1663  WRITE(*,*) 'sufficient function decrease'
1664 
1665  END IF
1666 
1667  EXIT optimal_step_search
1668 
1669  ELSE IF ( alam < alamin ) THEN
1670  ! convergence on Delta_x
1671 
1672  IF ( verbose_level .GE. 4 ) THEN
1673 
1674  WRITE(*,*) ' convergence on Delta_x',alam,alamin
1675 
1676  END IF
1677 
1678  qj_rel(:) = qj_rel_nr_old(:)
1679  scal_f = scal_f_old
1680  check = .true.
1681 
1682  EXIT optimal_step_search
1683 
1684  ! ELSE IF ( scal_f .LE. scal_f_old + ALF * alam * slope ) THEN
1685  ELSE
1686 
1687  IF ( alam .EQ. 1.d0 ) THEN
1688 
1689  tmplam = - slope / ( 2.0d0 * ( scal_f - scal_f_old - slope ) )
1690 
1691  ELSE
1692 
1693  rhs1 = scal_f - scal_f_old - alam*slope
1694  rhs2 = scal_f2 - scal_f_old - alam2*slope
1695 
1696  a = ( rhs1/alam**2.d0 - rhs2/alam2**2.d0 ) / ( alam - alam2 )
1697  b = ( -alam2*rhs1/alam**2 + alam*rhs2/alam2**2 ) / ( alam - alam2 )
1698 
1699  IF ( a .EQ. 0.d0 ) THEN
1700 
1701  tmplam = - slope / ( 2.0d0 * b )
1702 
1703  ELSE
1704 
1705  disc = b*b - 3.0d0*a*slope
1706 
1707  IF ( disc .LT. 0.d0 ) THEN
1708 
1709  tmplam = 0.5d0 * alam
1710 
1711  ELSE IF ( b .LE. 0.d0 ) THEN
1712 
1713  tmplam = ( - b + sqrt(disc) ) / ( 3.d0 * a )
1714 
1715  ELSE
1716 
1717  tmplam = - slope / ( b + sqrt(disc) )
1718 
1719  ENDIF
1720 
1721  END IF
1722 
1723  IF ( tmplam .GT. 0.5d0*alam ) tmplam = 0.5d0 * alam
1724 
1725  END IF
1726 
1727  END IF
1728 
1729  alam2 = alam
1730  scal_f2 = scal_f
1731  alam = max( tmplam , 0.5d0*alam )
1732 
1733  END DO optimal_step_search
1734 
1735  RETURN
1736 
1737  END SUBROUTINE lnsrch
1738 
1739  !******************************************************************************
1741  !
1755  !******************************************************************************
1756 
1757  SUBROUTINE eval_f( qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , f_nl , &
1758  scal_f )
1759 
1760  USE constitutive_2d, ONLY : eval_nonhyperbolic_terms
1761 
1762  IMPLICIT NONE
1763 
1764  REAL*8, INTENT(IN) :: qj(n_vars)
1765  REAL*8, INTENT(IN) :: qj_old(n_vars)
1766  REAL*8, INTENT(IN) :: a_tilde(n_rk)
1767  REAL*8, INTENT(IN) :: a_dirk(n_rk)
1768  REAL*8, INTENT(IN) :: a_diag
1769  REAL*8, INTENT(IN) :: coeff_f(n_eqns)
1770  REAL*8, INTENT(OUT) :: f_nl(n_eqns)
1771  REAL*8, INTENT(OUT) :: scal_f
1772 
1773  REAL*8 :: nh_term_impl(n_eqns)
1774  REAL*8 :: Rj(n_eqns)
1775 
1776  CALL eval_nonhyperbolic_terms( r_qj = qj , r_nh_term_impl=nh_term_impl )
1777 
1778  rj = ( matmul( divfluxj(1:n_eqns,1:i_rk-1)+expl_terms_j(1:n_eqns,1:i_rk-1), &
1779  a_tilde(1:i_rk-1) ) - matmul( nhj(1:n_eqns,1:i_rk-1) &
1780  + si_nhj(1:n_eqns,1:i_rk-1) , a_dirk(1:i_rk-1) ) ) &
1781  - a_diag * ( nh_term_impl + si_nhj(1:n_eqns,i_rk) )
1782 
1783  f_nl = qj - qj_old + dt * rj
1784 
1785  f_nl = coeff_f * f_nl
1786 
1787  scal_f = 0.5d0 * dot_product( f_nl , f_nl )
1788 
1789  RETURN
1790 
1791  END SUBROUTINE eval_f
1792 
1793  !******************************************************************************
1795  !
1798  !
1803  !
1807  !******************************************************************************
1808 
1809  SUBROUTINE eval_jacobian( qj_rel , qj_org , coeff_f, left_matrix)
1810 
1811  USE constitutive_2d, ONLY : eval_nonhyperbolic_terms
1812 
1813  IMPLICIT NONE
1814 
1815  REAL*8, INTENT(IN) :: qj_rel(n_vars)
1816  REAL*8, INTENT(IN) :: qj_org(n_vars)
1817  REAL*8, INTENT(IN) :: coeff_f(n_eqns)
1818 
1819  REAL*8, INTENT(OUT) :: left_matrix(n_eqns,n_vars)
1820 
1821  REAL*8 :: Jacob_relax(n_eqns,n_vars)
1822  COMPLEX*16 :: nh_terms_cmplx_impl(n_eqns)
1823  COMPLEX*16 :: qj_cmplx(n_vars) , qj_rel_cmplx(n_vars)
1824 
1825  REAL*8 :: h
1826 
1827  INTEGER :: i
1828 
1829  h = n_vars * epsilon(1.d0)
1830 
1831  ! initialize the matrix of the linearized system and the Jacobian
1832 
1833  left_matrix(1:n_eqns,1:n_vars) = 0.d0
1834  jacob_relax(1:n_eqns,1:n_vars) = 0.d0
1835 
1836  ! evaluate the jacobian of the non-hyperbolic terms
1837 
1838  DO i=1,n_vars
1839 
1840  left_matrix(i,i) = coeff_f(i) * qj_org(i)
1841 
1842  IF ( implicit_flag(i) ) THEN
1843 
1844  qj_rel_cmplx(1:n_vars) = qj_rel(1:n_vars)
1845  qj_rel_cmplx(i) = dcmplx(qj_rel(i), h)
1846 
1847  qj_cmplx = qj_rel_cmplx * qj_org
1848 
1849  CALL eval_nonhyperbolic_terms( c_qj = qj_cmplx , &
1850  c_nh_term_impl = nh_terms_cmplx_impl )
1851 
1852  jacob_relax(1:n_eqns,i) = coeff_f(i) * &
1853  aimag(nh_terms_cmplx_impl) / h
1854 
1855  left_matrix(1:n_eqns,i) = left_matrix(1:n_eqns,i) - dt * a_diag &
1856  * jacob_relax(1:n_eqns,i)
1857 
1858  END IF
1859 
1860  END DO
1861 
1862  RETURN
1863 
1864  END SUBROUTINE eval_jacobian
1865 
1866  !******************************************************************************
1868  !
1871  !
1873  !
1877  !******************************************************************************
1878 
1879  SUBROUTINE update_erosion_deposition_cell(dt)
1880 
1881  USE constitutive_2d, ONLY : erosion_coeff , settling_flag , rho_s
1882 
1883  USE geometry_2d, ONLY : deposit
1884 
1885  USE constitutive_2d, ONLY : eval_erosion_dep_term
1886  USE constitutive_2d, ONLY : eval_topo_term
1887 
1888  USE constitutive_2d, ONLY : qc_to_qp , mixt_var
1889  USE parameters_2d, ONLY : topo_change_flag
1890 
1891  IMPLICIT NONE
1892 
1893  REAL*8, INTENT(IN) :: dt
1894 
1895  REAL*8 :: erosion_term(n_solid)
1896  REAL*8 :: deposition_term(n_solid)
1897  REAL*8 :: eqns_term(n_eqns)
1898  REAL*8 :: deposit_term(n_solid)
1899 
1900  REAL*8 :: r_Ri , r_rho_m
1901  REAL*8 :: r_rho_c
1902  REAL*8 :: r_red_grav
1903 
1904  INTEGER :: j ,k , l
1905 
1906  IF ( ( sum(erosion_coeff) .EQ. 0.d0 ) .AND. ( .NOT. settling_flag ) ) RETURN
1907 
1908  DO l = 1,solve_cells
1909 
1910  j = j_cent(l)
1911  k = k_cent(l)
1912 
1913  CALL qc_to_qp(q(1:n_vars,j,k) , qp(1:n_vars+2,j,k) )
1914 
1915  CALL eval_erosion_dep_term( qp(1:n_vars+2,j,k) , dt , &
1916  erosion_term(1:n_solid) , deposition_term(1:n_solid) )
1917 
1918  ! Limit the deposition during a single time step
1919  deposition_term(1:n_solid) = max(0.d0,min( deposition_term(1:n_solid), &
1920  q(5:4+n_solid,j,k) / ( rho_s(1:n_solid) * dt ) ))
1921 
1922  ! Compute the source terms for the equations
1923  CALL eval_topo_term( qp(1:n_vars+2,j,k) , deposition_term , &
1924  erosion_term , eqns_term , deposit_term )
1925 
1926  IF ( verbose_level .GE. 2 ) THEN
1927 
1928  WRITE(*,*) 'before update erosion/deposition: j,k,q(:,j,k),B(j,k)', &
1929  j,k,q(:,j,k),b_cent(j,k)
1930 
1931  END IF
1932 
1933  ! Update the solution with erosion/deposition terms
1934  q(1:n_eqns,j,k) = q(1:n_eqns,j,k) + dt * eqns_term(1:n_eqns)
1935 
1936  deposit(j,k,1:n_solid) = deposit(j,k,1:n_solid)+dt*deposit_term(1:n_solid)
1937 
1938  ! Update the topography with erosion/deposition terms
1939  IF ( topo_change_flag ) THEN
1940 
1941  b_cent(j,k) = b_cent(j,k) + dt * sum( deposit_term(1:n_solid) )
1942 
1943  END IF
1944 
1945  ! Check for negative thickness
1946  IF ( q(1,j,k) .LT. 0.d0 ) THEN
1947 
1948  IF ( q(1,j,k) .GT. -1.d-10 ) THEN
1949 
1950  q(1:n_vars,j,k) = 0.d0
1951 
1952  ELSE
1953 
1954  WRITE(*,*) 'j,k',j,k
1955  WRITE(*,*) 'dt',dt
1956  WRITE(*,*) 'before erosion'
1957  WRITE(*,*) 'qp',qp(1:n_eqns+2,j,k)
1958  WRITE(*,*) 'deposition_term',deposition_term
1959  WRITE(*,*) 'erosion_term',erosion_term
1960  CALL qc_to_qp(q(1:n_vars,j,k) , qp(1:n_vars+2,j,k))
1961  WRITE(*,*) 'qp',qp(1:n_eqns+2,j,k)
1962  READ(*,*)
1963 
1964  END IF
1965 
1966  END IF
1967 
1968  CALL qc_to_qp(q(1:n_vars,j,k) , qp(1:n_vars+2,j,k) )
1969  CALL mixt_var(qp(1:n_vars+2,j,k),r_ri,r_rho_m,r_rho_c,r_red_grav)
1970 
1971  IF ( r_red_grav .LE. 0.d0 ) THEN
1972 
1973  q(1:n_vars,j,k) = 0.d0
1974 
1975  END IF
1976 
1977  END DO
1978 
1979  RETURN
1980 
1981  END SUBROUTINE update_erosion_deposition_cell
1982 
1983  !******************************************************************************
1985  !
1988  !
1992  !
1996  !******************************************************************************
1997 
1998  SUBROUTINE eval_explicit_terms( qc_expl , qp_expl , expl_terms )
1999 
2000  USE constitutive_2d, ONLY : eval_expl_terms
2001 
2002  IMPLICIT NONE
2003 
2004  REAL*8, INTENT(IN) :: qc_expl(n_vars,comp_cells_x,comp_cells_y)
2005  REAL*8, INTENT(IN) :: qp_expl(n_vars+2,comp_cells_x,comp_cells_y)
2006  REAL*8, INTENT(OUT) :: expl_terms(n_eqns,comp_cells_x,comp_cells_y)
2007 
2008  REAL*8 :: qcj(n_vars)
2009  REAL*8 :: qpj(n_vars+2)
2010  REAL*8 :: expl_forces_term(n_eqns)
2011 
2012  INTEGER :: j,k,l
2013 
2014  expl_terms = 0.d0
2015 
2016  DO l = 1,solve_cells
2017 
2018  j = j_cent(l)
2019  k = k_cent(l)
2020 
2021  qcj(1:n_vars) = qc_expl(1:n_vars,j,k)
2022  qpj(1:n_vars+2) = qp_expl(1:n_vars+2,j,k)
2023 
2024  CALL eval_expl_terms( b_prime_x(j,k), b_prime_y(j,k), source_xy(j,k) ,&
2025  qpj(1:n_vars+2) , qcj(1:n_vars) , expl_forces_term )
2026 
2027  expl_terms(1:n_eqns,j,k) = expl_forces_term
2028 
2029  END DO
2030 
2031  RETURN
2032 
2033  END SUBROUTINE eval_explicit_terms
2034 
2035  !******************************************************************************
2037  !
2042  !
2046  !
2050  !******************************************************************************
2051 
2052  SUBROUTINE eval_hyperbolic_terms( q_expl , qp_expl , divFlux )
2053 
2054  ! External variables
2055  USE geometry_2d, ONLY : dx,dy
2056  USE parameters_2d, ONLY : solver_scheme
2057 
2058  IMPLICIT NONE
2059 
2060  REAL*8, INTENT(IN) :: q_expl(n_vars,comp_cells_x,comp_cells_y)
2061  REAL*8, INTENT(IN) :: qp_expl(n_vars+2,comp_cells_x,comp_cells_y)
2062  REAL*8, INTENT(OUT) :: divFlux(n_eqns,comp_cells_x,comp_cells_y)
2063 
2064  REAL*8 :: h_new , h_old
2065 
2066  ! REAL*8 :: tcpu0,tcpu1,tcpu2,tcpu3,tcpu4
2067 
2068  INTEGER :: l , i, j, k
2069 
2070  ! CALL cpu_time(tcpu0)
2071 
2072  ! Linear reconstruction of the physical variables at the interfaces
2073  CALL reconstruction(q_expl,qp_expl)
2074 
2075  ! CALL cpu_time(tcpu1)
2076  !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu1-tcpu0,'seconds'
2077 
2078  ! Evaluation of the maximum local speeds at the interfaces
2079  CALL eval_speeds
2080 
2081  ! CALL cpu_time(tcpu2)
2082  !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu2-tcpu1,'seconds'
2083 
2084  ! Evaluation of the numerical fluxes
2085  SELECT CASE ( solver_scheme )
2086 
2087  CASE ("LxF")
2088 
2089  CALL eval_flux_lxf
2090 
2091  CASE ("GFORCE")
2092 
2093  CALL eval_flux_gforce
2094 
2095  CASE ("KT")
2096 
2097  CALL eval_flux_kt
2098 
2099  CASE ("UP")
2100 
2101  CALL eval_flux_up
2102 
2103  END SELECT
2104 
2105  ! CALL cpu_time(tcpu3)
2106  !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu3-tcpu2,'seconds'
2107 
2108  ! Advance in time the solution
2109  DO l = 1,solve_cells
2110 
2111  j = j_cent(l)
2112  k = k_cent(l)
2113 
2114  DO i=1,n_eqns
2115 
2116  divflux(i,j,k) = 0.d0
2117 
2118  IF ( comp_cells_x .GT. 1 ) THEN
2119 
2120  divflux(i,j,k) = divflux(i,j,k) + &
2121  ( h_interface_x(i,j+1,k) - h_interface_x(i,j,k) ) / dx
2122 
2123  END IF
2124 
2125  IF ( comp_cells_y .GT. 1 ) THEN
2126 
2127  divflux(i,j,k) = divflux(i,j,k) + &
2128  ( h_interface_y(i,j,k+1) - h_interface_y(i,j,k) ) / dy
2129 
2130  END IF
2131 
2132  END DO
2133 
2134  h_old = q_expl(1,j,k)
2135  h_new = h_old - dt * divflux(1,j,k)
2136 
2137  END DO
2138 
2139  ! CALL cpu_time(tcpu4)
2140  !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu4-tcpu4,'seconds'
2141 
2142  RETURN
2143 
2144  END SUBROUTINE eval_hyperbolic_terms
2145 
2146  !******************************************************************************
2148  !
2154  !******************************************************************************
2155 
2156  SUBROUTINE eval_flux_up
2157 
2158  ! External procedures
2159  USE constitutive_2d, ONLY : eval_fluxes
2160 
2161  IMPLICIT NONE
2162 
2163  REAL*8 :: fluxL(n_eqns)
2164  REAL*8 :: fluxR(n_eqns)
2165  REAL*8 :: fluxB(n_eqns)
2166  REAL*8 :: fluxT(n_eqns)
2167 
2168  INTEGER :: j,k,l
2169 
2170  h_interface_x = 0.d0
2171  h_interface_y = 0.d0
2172 
2173  IF ( comp_cells_x .GT. 1 ) THEN
2174 
2175  DO l = 1,solve_interfaces_x
2176 
2177  j = j_stag_x(l)
2178  k = k_stag_x(l)
2179 
2180 
2181  CALL eval_fluxes( q_interfacel(1:n_vars,j,k) , &
2182  qp_interfacel(1:n_vars+2,j,k) , 1 , fluxl)
2183 
2184  CALL eval_fluxes( q_interfacer(1:n_vars,j,k) , &
2185  qp_interfacer(1:n_vars+2,j,k) , 1 , fluxr)
2186 
2187  IF ( ( q_interfacel(2,j,k) .GT. 0.d0 ) .AND. &
2188  ( q_interfacer(2,j,k) .GE. 0.d0 ) ) THEN
2189 
2190  h_interface_x(:,j,k) = fluxl
2191 
2192  ELSEIF ( ( q_interfacel(2,j,k) .LE. 0.d0 ) .AND. &
2193  ( q_interfacer(2,j,k) .LT. 0.d0 ) ) THEN
2194 
2195  h_interface_x(:,j,k) = fluxr
2196 
2197  ELSE
2198 
2199  h_interface_x(:,j,k) = 0.5d0 * ( fluxl + fluxr )
2200 
2201  END IF
2202 
2203  IF ( ( q_interfacel(2,j,k) .EQ. 0.d0 ) .AND. &
2204  ( q_interfacer(2,j,k) .EQ. 0.d0 ) ) THEN
2205 
2206  h_interface_x(1,j,k) = 0.d0
2207  h_interface_x(4:n_vars,j,k) = 0.d0
2208 
2209  END IF
2210 
2211  END DO
2212 
2213  END IF
2214 
2215 
2216  IF ( comp_cells_y .GT. 1 ) THEN
2217 
2218  DO l = 1,solve_interfaces_y
2219 
2220  j = j_stag_y(l)
2221  k = k_stag_y(l)
2222 
2223  CALL eval_fluxes( q_interfaceb(1:n_vars,j,k) , &
2224  qp_interfaceb(1:n_vars+2,j,k) , 2 , fluxb)
2225 
2226  CALL eval_fluxes( q_interfacet(1:n_vars,j,k) , &
2227  qp_interfacet(1:n_vars+2,j,k) ,2 , fluxt)
2228 
2229  IF ( ( q_interfaceb(3,j,k) .GT. 0.d0 ) .AND. &
2230  ( q_interfacet(3,j,k) .GE. 0.d0 ) ) THEN
2231 
2232  h_interface_y(:,j,k) = fluxb
2233 
2234  ELSEIF ( ( q_interfaceb(3,j,k) .LE. 0.d0 ) .AND. &
2235  ( q_interfacet(3,j,k) .LT. 0.d0 ) ) THEN
2236 
2237  h_interface_y(:,j,k) = fluxt
2238 
2239  ELSE
2240 
2241  h_interface_y(:,j,k) = 0.5d0 * ( fluxb + fluxt )
2242 
2243  END IF
2244 
2245  ! In the equation for mass and for trasnport (T,alphas) if the
2246  ! velocities at the interfaces are null, then the flux is null
2247  IF ( ( q_interfaceb(3,j,k) .EQ. 0.d0 ) .AND. &
2248  ( q_interfacet(3,j,k) .EQ. 0.d0 ) ) THEN
2249 
2250  h_interface_y(1,j,k) = 0.d0
2251  h_interface_y(4:n_vars,j,k) = 0.d0
2252 
2253  END IF
2254 
2255  END DO
2256 
2257  END IF
2258 
2259  RETURN
2260 
2261  END SUBROUTINE eval_flux_up
2262 
2263 
2264  !******************************************************************************
2266  !
2272  !******************************************************************************
2273 
2274  SUBROUTINE eval_flux_kt
2275 
2276  ! External procedures
2277  USE constitutive_2d, ONLY : eval_fluxes
2278 
2279  IMPLICIT NONE
2280 
2281  REAL*8 :: fluxL(n_eqns)
2282  REAL*8 :: fluxR(n_eqns)
2283  REAL*8 :: fluxB(n_eqns)
2284  REAL*8 :: fluxT(n_eqns)
2285 
2286  REAL*8 :: flux_avg_x(n_eqns)
2287  REAL*8 :: flux_avg_y(n_eqns)
2288 
2289  INTEGER :: i,j,k,l
2290 
2291  h_interface_x = 0.d0
2292  h_interface_y = 0.d0
2293 
2294  IF ( comp_cells_x .GT. 1 ) THEN
2295 
2296  DO l = 1,solve_interfaces_x
2297 
2298  j = j_stag_x(l)
2299  k = k_stag_x(l)
2300 
2301 
2302  CALL eval_fluxes( q_interfacel(1:n_vars,j,k) , &
2303  qp_interfacel(1:n_vars+2,j,k) , 1 , fluxl)
2304 
2305  CALL eval_fluxes( q_interfacer(1:n_vars,j,k) , &
2306  qp_interfacer(1:n_vars+2,j,k) , 1 , fluxr)
2307 
2308  CALL average_kt( a_interface_xneg(:,j,k), a_interface_xpos(:,j,k) , &
2309  fluxl , fluxr , flux_avg_x )
2310 
2311  eqns_loop:DO i=1,n_eqns
2312 
2313  IF ( a_interface_xneg(i,j,k) .EQ. a_interface_xpos(i,j,k) ) THEN
2314 
2315  h_interface_x(i,j,k) = 0.d0
2316 
2317  ELSE
2318 
2319  h_interface_x(i,j,k) = flux_avg_x(i) &
2320  + ( a_interface_xpos(i,j,k) * a_interface_xneg(i,j,k) ) &
2321  / ( a_interface_xpos(i,j,k) - a_interface_xneg(i,j,k) ) &
2322  * ( q_interfacer(i,j,k) - q_interfacel(i,j,k) )
2323 
2324  END IF
2325 
2326  ENDDO eqns_loop
2327 
2328  ! In the equation for mass and for trasnport (T,alphas) if the
2329  ! velocities at the interfaces are null, then the flux is null
2330  IF ( ( qp_interfacel(2,j,k) .EQ. 0.d0 ) .AND. &
2331  ( qp_interfacer(2,j,k) .EQ. 0.d0 ) ) THEN
2332 
2333  h_interface_x(1,j,k) = 0.d0
2334  h_interface_x(4:n_vars,j,k) = 0.d0
2335 
2336  END IF
2337 
2338  END DO
2339 
2340  END IF
2341 
2342 
2343  IF ( comp_cells_y .GT. 1 ) THEN
2344 
2345  DO l = 1,solve_interfaces_y
2346 
2347  j = j_stag_y(l)
2348  k = k_stag_y(l)
2349 
2350  CALL eval_fluxes( q_interfaceb(1:n_vars,j,k) , &
2351  qp_interfaceb(1:n_vars+2,j,k) , 2 , fluxb)
2352 
2353  CALL eval_fluxes( q_interfacet(1:n_vars,j,k) , &
2354  qp_interfacet(1:n_vars+2,j,k) , 2 , fluxt)
2355 
2356  CALL average_kt( a_interface_yneg(:,j,k) , &
2357  a_interface_ypos(:,j,k) , fluxb , fluxt , flux_avg_y )
2358 
2359  DO i=1,n_eqns
2360 
2361  IF ( a_interface_yneg(i,j,k) .EQ. a_interface_ypos(i,j,k) ) THEN
2362 
2363  h_interface_y(i,j,k) = 0.d0
2364 
2365  ELSE
2366 
2367  h_interface_y(i,j,k) = flux_avg_y(i) &
2368  + ( a_interface_ypos(i,j,k) * a_interface_yneg(i,j,k) ) &
2369  / ( a_interface_ypos(i,j,k) - a_interface_yneg(i,j,k) ) &
2370  * ( q_interfacet(i,j,k) - q_interfaceb(i,j,k) )
2371 
2372  END IF
2373 
2374  END DO
2375 
2376  ! In the equation for mass and for trasnport (T,alphas) if the
2377  ! velocities at the interfaces are null, then the flux is null
2378  IF ( ( q_interfaceb(3,j,k) .EQ. 0.d0 ) .AND. &
2379  ( q_interfacet(3,j,k) .EQ. 0.d0 ) ) THEN
2380 
2381  h_interface_y(1,j,k) = 0.d0
2382  h_interface_y(4:n_vars,j,k) = 0.d0
2383 
2384  END IF
2385 
2386  END DO
2387 
2388  END IF
2389 
2390  RETURN
2391 
2392  END SUBROUTINE eval_flux_kt
2393 
2394  !******************************************************************************
2396  !
2407  !******************************************************************************
2408 
2409  SUBROUTINE average_kt( a1 , a2 , w1 , w2 , w_avg )
2410 
2411  IMPLICIT NONE
2412 
2413  REAL*8, INTENT(IN) :: a1(:) , a2(:)
2414  REAL*8, INTENT(IN) :: w1(:) , w2(:)
2415  REAL*8, INTENT(OUT) :: w_avg(:)
2416 
2417  INTEGER :: n
2418  INTEGER :: i
2419 
2420  n = SIZE( a1 )
2421 
2422  DO i=1,n
2423 
2424  IF ( a1(i) .EQ. a2(i) ) THEN
2425 
2426  w_avg(i) = 0.5d0 * ( w1(i) + w2(i) )
2427  w_avg(i) = 0.d0
2428 
2429  ELSE
2430 
2431  w_avg(i) = ( a2(i) * w1(i) - a1(i) * w2(i) ) / ( a2(i) - a1(i) )
2432 
2433  END IF
2434 
2435  END DO
2436 
2437  RETURN
2438 
2439  END SUBROUTINE average_kt
2440 
2441  !******************************************************************************
2446  !******************************************************************************
2447 
2448  SUBROUTINE eval_flux_gforce
2449 
2450  ! to be implemented
2451  WRITE(*,*) 'method not yet implemented in 2-d case'
2452 
2453  END SUBROUTINE eval_flux_gforce
2454 
2455  !******************************************************************************
2460  !******************************************************************************
2461 
2462  SUBROUTINE eval_flux_lxf
2463 
2464  ! to be implemented
2465  WRITE(*,*) 'method not yet implemented in 2-d case'
2466 
2467  END SUBROUTINE eval_flux_lxf
2468 
2469 
2470  !******************************************************************************
2472  !
2483  !******************************************************************************
2484 
2485  SUBROUTINE reconstruction(q_expl,qp_expl)
2486 
2487  ! External procedures
2488  USE constitutive_2d, ONLY : qc_to_qp , qp_to_qc ,qp_to_qp2
2489  USE constitutive_2d, ONLY : eval_source_bdry
2490  USE parameters_2d, ONLY : limiter
2491 
2492  ! External variables
2493  USE geometry_2d, ONLY : x_comp , x_stag , y_comp , y_stag , dx , dx2 , dy , &
2494  dy2
2495 
2496  USE geometry_2d, ONLY : b_cent
2497  USE geometry_2d, ONLY : b_interfacel , b_interfacer
2498  USE geometry_2d, ONLY : b_interfacet , b_interfaceb
2499 
2500  USE geometry_2d, ONLY : sourcew , sourcee , sourcen , sources
2501  USE geometry_2d, ONLY : sourcew_vect_x , sourcew_vect_y
2502  USE geometry_2d, ONLY : sourcee_vect_x , sourcee_vect_y
2503  USE geometry_2d, ONLY : sourcen_vect_x , sourcen_vect_y
2504  USE geometry_2d, ONLY : sources_vect_x , sources_vect_y
2505 
2506  USE parameters_2d, ONLY : reconstr_coeff
2507 
2508  USE geometry_2d, ONLY : minmod
2509 
2510  IMPLICIT NONE
2511 
2512  REAL*8, INTENT(IN) :: q_expl(:,:,:)
2513  REAL*8, INTENT(IN) :: qp_expl(:,:,:)
2514 
2515  REAL*8, ALLOCATABLE :: qrec(:,:,:)
2516  REAL*8, ALLOCATABLE :: qrecW(:)
2517  REAL*8, ALLOCATABLE :: qrecE(:)
2518  REAL*8, ALLOCATABLE :: qrecS(:)
2519  REAL*8, ALLOCATABLE :: qrecN(:)
2520 
2521  REAL*8, ALLOCATABLE :: source_bdry(:)
2522  REAL*8, ALLOCATABLE :: qrec_prime_x(:)
2523  REAL*8, ALLOCATABLE :: qrec_prime_y(:)
2524 
2525  REAL*8 :: qp2recW(3) , qp2recE(3)
2526  REAL*8 :: qp2recS(3) , qp2recN(3)
2527 
2528  REAL*8 :: qrec_stencil(3)
2529  REAL*8 :: x_stencil(3)
2530  REAL*8 :: y_stencil(3)
2531 
2532  INTEGER :: l,j,k
2533  INTEGER :: i
2534 
2535 
2536  ALLOCATE ( qrec(n_vars+2,comp_cells_x,comp_cells_y) )
2537  ALLOCATE ( qrecw(n_vars+2) )
2538  ALLOCATE ( qrece(n_vars+2) )
2539  ALLOCATE ( qrecs(n_vars+2) )
2540  ALLOCATE ( qrecn(n_vars+2) )
2541 
2542  ALLOCATE ( source_bdry(n_vars+2) )
2543 
2544  ALLOCATE ( qrec_prime_x(n_vars+2) )
2545  ALLOCATE ( qrec_prime_y(n_vars+2) )
2546 
2547  ! Compute the variable to reconstruct (phys or cons)
2548  DO k = 1,comp_cells_y
2549 
2550  DO j = 1,comp_cells_x
2551 
2552  qrec(1:n_vars+2,j,k) = qp_expl(1:n_vars+2,j,k)
2553 
2554  IF ( sum(qrec(5:4+n_solid,j,k)) .GT. 1.d0 ) THEN
2555 
2556  WRITE(*,*) 'reconstruction: j,k',j,k
2557  WRITE(*,*) 'qrec(5:n_solid,j,k)',qrec(5:4+n_solid,j,k)
2558  WRITE(*,*) 'q(1:n_vars,j,k)',q_expl(1:n_vars,j,k)
2559  WRITE(*,*) 'B_cent(j,k)', b_cent(j,k)
2560  READ(*,*)
2561 
2562  END IF
2563 
2564  END DO
2565 
2566  END DO
2567 
2568  ! Linear reconstruction
2569  DO l = 1,solve_cells
2570 
2571  j = j_cent(l)
2572  k = k_cent(l)
2573 
2574 
2575  qrecw(1:n_vars+2) = qp_expl(1:n_vars+2,j,k)
2576  qrece(1:n_vars+2) = qp_expl(1:n_vars+2,j,k)
2577  qrecs(1:n_vars+2) = qp_expl(1:n_vars+2,j,k)
2578  qrecn(1:n_vars+2) = qp_expl(1:n_vars+2,j,k)
2579 
2580  vars_loop:DO i=1,n_vars
2581 
2582  ! x direction
2583  check_comp_cells_x:IF ( comp_cells_x .GT. 1 ) THEN
2584 
2585  ! west boundary
2586  check_x_boundary:IF (j.EQ.1) THEN
2587 
2588  IF ( bcw(i)%flag .EQ. 0 ) THEN
2589 
2590  x_stencil(1) = x_stag(1)
2591  x_stencil(2:3) = x_comp(1:2)
2592 
2593  qrec_stencil(1) = bcw(i)%value
2594  qrec_stencil(2:3) = qrec(i,1:2,k)
2595 
2596  CALL limit( qrec_stencil , x_stencil , limiter(i) , &
2597  qrec_prime_x(i) )
2598 
2599  ELSEIF ( bcw(i)%flag .EQ. 1 ) THEN
2600 
2601  qrec_prime_x(i) = bcw(i)%value
2602 
2603  ELSEIF ( bcw(i)%flag .EQ. 2 ) THEN
2604 
2605  qrec_prime_x(i) = ( qrec(i,2,k) - qrec(i,1,k) ) / dx
2606 
2607  END IF
2608 
2609  !east boundary
2610  ELSEIF (j.EQ.comp_cells_x) THEN
2611 
2612  IF ( bce(i)%flag .EQ. 0 ) THEN
2613 
2614  qrec_stencil(3) = bce(i)%value
2615  qrec_stencil(1:2)= qrec(i,comp_cells_x-1:comp_cells_x,k)
2616 
2617  x_stencil(3) = x_stag(comp_interfaces_x)
2618  x_stencil(1:2) = x_comp(comp_cells_x-1:comp_cells_x)
2619 
2620  CALL limit( qrec_stencil , x_stencil , limiter(i) , &
2621  qrec_prime_x(i) )
2622 
2623  ELSEIF ( bce(i)%flag .EQ. 1 ) THEN
2624 
2625  qrec_prime_x(i) = bce(i)%value
2626 
2627  ELSEIF ( bce(i)%flag .EQ. 2 ) THEN
2628 
2629  qrec_prime_x(i) = ( qrec(i,comp_cells_x,k) - &
2630  qrec(i,comp_cells_x-1,k) ) / dx
2631 
2632  END IF
2633 
2634  ! internal x cells
2635  ELSE
2636 
2637  x_stencil(1:3) = x_comp(j-1:j+1)
2638  qrec_stencil(1:3) = qrec(i,j-1:j+1,k)
2639 
2640  ! correction for radial source inlet x-interfaces values
2641  ! used for the linear reconstruction
2642  IF ( radial_source_flag .AND. ( source_cell(j,k).EQ.2 ) ) THEN
2643 
2644  IF ( sourcee(j,k) ) THEN
2645 
2646  CALL eval_source_bdry( t, sourcee_vect_x(j,k) , &
2647  sourcee_vect_y(j,k) , source_bdry )
2648 
2649  x_stencil(3) = x_stag(j+1)
2650  qrec_stencil(3) = source_bdry(i)
2651 
2652  ELSEIF ( sourcew(j,k) ) THEN
2653 
2654  CALL eval_source_bdry( t , sourcew_vect_x(j,k) , &
2655  sourcew_vect_y(j,k) , source_bdry )
2656 
2657  x_stencil(1) = x_stag(j)
2658  qrec_stencil(1) = source_bdry(i)
2659 
2660  END IF
2661 
2662  END IF
2663 
2664  CALL limit( qrec_stencil , x_stencil , limiter(i) , &
2665  qrec_prime_x(i) )
2666 
2667  ENDIF check_x_boundary
2668 
2669  qrecw(i) = qrec(i,j,k) - reconstr_coeff* dx2 * qrec_prime_x(i)
2670  qrece(i) = qrec(i,j,k) + reconstr_coeff* dx2 * qrec_prime_x(i)
2671 
2672  END IF check_comp_cells_x
2673 
2674  ! y-direction
2675  check_comp_cells_y:IF ( comp_cells_y .GT. 1 ) THEN
2676 
2677  ! South boundary
2678  check_y_boundary:IF (k.EQ.1) THEN
2679 
2680  IF ( bcs(i)%flag .EQ. 0 ) THEN
2681 
2682  qrec_stencil(1) = bcs(i)%value
2683  qrec_stencil(2:3) = qrec(i,j,1:2)
2684 
2685  y_stencil(1) = y_stag(1)
2686  y_stencil(2:3) = y_comp(1:2)
2687 
2688  CALL limit( qrec_stencil , y_stencil , limiter(i) , &
2689  qrec_prime_y(i) )
2690 
2691  ELSEIF ( bcs(i)%flag .EQ. 1 ) THEN
2692 
2693  qrec_prime_y(i) = bcs(i)%value
2694 
2695  ELSEIF ( bcs(i)%flag .EQ. 2 ) THEN
2696 
2697  qrec_prime_y(i) = ( qrec(i,j,2) - qrec(i,j,1) ) / dy
2698 
2699  END IF
2700 
2701  ! North boundary
2702  ELSEIF ( k .EQ. comp_cells_y ) THEN
2703 
2704  IF ( bcn(i)%flag .EQ. 0 ) THEN
2705 
2706  qrec_stencil(3) = bcn(i)%value
2707  qrec_stencil(1:2)= qrec(i,j,comp_cells_y-1:comp_cells_y)
2708 
2709  y_stencil(3) = y_stag(comp_interfaces_y)
2710  y_stencil(1:2) = y_comp(comp_cells_y-1:comp_cells_y)
2711 
2712  CALL limit( qrec_stencil , y_stencil , limiter(i) , &
2713  qrec_prime_y(i) )
2714 
2715  ELSEIF ( bcn(i)%flag .EQ. 1 ) THEN
2716 
2717  qrec_prime_y(i) = bcn(i)%value
2718 
2719  ELSEIF ( bcn(i)%flag .EQ. 2 ) THEN
2720 
2721  qrec_prime_y(i) = ( qrec(i,j,comp_cells_y) - &
2722  qrec(i,j,comp_cells_y-1) ) / dy
2723 
2724  END IF
2725 
2726  ! Internal y cells
2727  ELSE
2728 
2729  y_stencil(1:3) = y_comp(k-1:k+1)
2730  qrec_stencil = qrec(i,j,k-1:k+1)
2731 
2732  ! correction for radial source inlet y-interfaces
2733  ! used for the linear reconstruction
2734  IF ( radial_source_flag .AND. ( source_cell(j,k).EQ.2 ) ) THEN
2735 
2736  IF ( sources(j,k) ) THEN
2737 
2738  CALL eval_source_bdry( t, sources_vect_x(j,k) , &
2739  sources_vect_y(j,k) , source_bdry )
2740 
2741  y_stencil(1) = y_stag(k)
2742  qrec_stencil(1) = source_bdry(i)
2743 
2744  ELSEIF ( sourcen(j,k) ) THEN
2745 
2746  CALL eval_source_bdry( t, sourcen_vect_x(j,k) , &
2747  sourcen_vect_y(j,k) , source_bdry )
2748 
2749  x_stencil(3) = y_stag(k+1)
2750  qrec_stencil(3) = source_bdry(i)
2751 
2752  END IF
2753 
2754  END IF
2755 
2756  CALL limit( qrec_stencil , y_stencil , limiter(i) , &
2757  qrec_prime_y(i) )
2758 
2759  ENDIF check_y_boundary
2760 
2761  qrecs(i) = qrec(i,j,k) - reconstr_coeff* dy2 * qrec_prime_y(i)
2762  qrecn(i) = qrec(i,j,k) + reconstr_coeff* dy2 * qrec_prime_y(i)
2763 
2764  ENDIF check_comp_cells_y
2765 
2766  ENDDO vars_loop
2767 
2768  add_vars_loop:DO i=n_vars+1,n_vars+2
2769  ! reconstruction on u and v with same limiters of hu,hv
2770 
2771  ! x direction
2772  check_comp_cells_x2:IF ( comp_cells_x .GT. 1 ) THEN
2773 
2774  IF ( (j .GT. 1) .AND. (j .LT. comp_cells_x) ) THEN
2775 
2776  x_stencil(1:3) = x_comp(j-1:j+1)
2777  qrec_stencil(1:3) = qrec(i,j-1:j+1,k)
2778 
2779  ! correction for radial source inlet x-interfaces values
2780  ! used for the linear reconstruction
2781  IF ( radial_source_flag .AND. ( source_cell(j,k).EQ.2 ) ) THEN
2782 
2783  IF ( sourcee(j,k) ) THEN
2784 
2785  CALL eval_source_bdry( t, sourcee_vect_x(j,k) , &
2786  sourcee_vect_y(j,k) , source_bdry )
2787 
2788  x_stencil(3) = x_stag(j+1)
2789  qrec_stencil(3) = source_bdry(i)
2790 
2791  ELSEIF ( sourcew(j,k) ) THEN
2792 
2793  CALL eval_source_bdry( t , sourcew_vect_x(j,k) , &
2794  sourcew_vect_y(j,k) , source_bdry )
2795 
2796  x_stencil(1) = x_stag(j)
2797  qrec_stencil(1) = source_bdry(i)
2798 
2799  END IF
2800 
2801  END IF
2802 
2803  CALL limit( qrec_stencil , x_stencil , limiter(i) , &
2804  qrec_prime_x(i) )
2805 
2806  qrecw(i) = qrec(i,j,k) - reconstr_coeff*dx2*qrec_prime_x(i)
2807  qrece(i) = qrec(i,j,k) + reconstr_coeff*dx2*qrec_prime_x(i)
2808 
2809  END IF
2810 
2811  END IF check_comp_cells_x2
2812 
2813  ! y-direction
2814  check_comp_cells_y2:IF ( comp_cells_y .GT. 1 ) THEN
2815 
2816  IF ( ( k .GT. 1 ) .AND. ( k .LT. comp_cells_y ) ) THEN
2817 
2818  y_stencil(1:3) = y_comp(k-1:k+1)
2819  qrec_stencil = qrec(i,j,k-1:k+1)
2820 
2821  ! correction for radial source inlet y-interfaces
2822  ! used for the linear reconstruction
2823  IF ( radial_source_flag .AND. ( source_cell(j,k).EQ.2 ) ) THEN
2824 
2825  IF ( sources(j,k) ) THEN
2826 
2827  CALL eval_source_bdry( t, sources_vect_x(j,k) , &
2828  sources_vect_y(j,k) , source_bdry )
2829 
2830  y_stencil(1) = y_stag(k)
2831  qrec_stencil(1) = source_bdry(i)
2832 
2833  ELSEIF ( sourcen(j,k) ) THEN
2834 
2835  CALL eval_source_bdry( t, sourcen_vect_x(j,k) , &
2836  sourcen_vect_y(j,k) , source_bdry )
2837 
2838  x_stencil(3) = y_stag(k+1)
2839  qrec_stencil(3) = source_bdry(i)
2840 
2841  END IF
2842 
2843  END IF
2844 
2845  CALL limit( qrec_stencil , y_stencil , limiter(i) , &
2846  qrec_prime_y(i) )
2847 
2848  qrecs(i) = qrec(i,j,k) - reconstr_coeff*dy2*qrec_prime_y(i)
2849  qrecn(i) = qrec(i,j,k) + reconstr_coeff*dy2*qrec_prime_y(i)
2850 
2851  ENDIF
2852 
2853  ENDIF check_comp_cells_y2
2854 
2855  ENDDO add_vars_loop
2856 
2857 
2858  IF ( comp_cells_x .GT. 1 ) THEN
2859 
2860  IF ( ( j .GT. 1 ) .AND. ( j .LT. comp_cells_x ) ) THEN
2861 
2862  IF ( q_expl(1,j,k) .EQ. 0.d0 ) THEN
2863 
2864  IF ( ( .NOT. radial_source_flag ) .OR. &
2865  ( ( radial_source_flag ) .AND. &
2866  ( source_cell(j,k) .EQ. 0 ) ) ) THEN
2867 
2868  ! In the internal cell, if thickness h is 0 at the center
2869  ! of the cell, then all the variables are 0 at the center
2870  ! and at the interfaces (no conversion back is needed from
2871  ! reconstructed to conservative)
2872  q_interfacer(:,j,k) = 0.d0
2873  q_interfacel(:,j+1,k) = 0.d0
2874 
2875  qp_interfacer(1:3,j,k) = 0.d0
2876  qp_interfacer(4:n_vars,j,k) = qrecw(4:n_vars)
2877  qp_interfacer(n_vars+1:n_vars+2,j,k) = 0.d0
2878 
2879  qp_interfacel(1:3,j+1,k) = 0.d0
2880  qp_interfacel(4:n_vars,j+1,k) = qrece(4:n_vars)
2881  qp_interfacel(n_vars+1:n_vars+2,j+1,k) = 0.d0
2882 
2883  END IF
2884 
2885  END IF
2886 
2887  END IF
2888 
2889  IF ( j.EQ.1 ) THEN
2890 
2891  ! Dirichelet boundary condition at the west of the domain
2892  DO i=1,n_vars
2893 
2894  IF ( bcw(i)%flag .EQ. 0 ) THEN
2895 
2896  qrecw(i) = bcw(i)%value
2897 
2898  END IF
2899 
2900  ENDDO
2901 
2902  DO i=n_vars+1,n_vars+2
2903 
2904  CALL qp_to_qp2( qrecw(1:n_vars) , b_cent(j,k) , qp2recw )
2905  qrecw(i) = qp2recw(i-n_vars+1)
2906 
2907  CALL qp_to_qp2( qrece(1:n_vars) , b_cent(j,k) , qp2rece )
2908  qrece(i) = qp2rece(i-n_vars+1)
2909 
2910  END DO
2911 
2912  ELSEIF ( j.EQ.comp_cells_x ) THEN
2913 
2914  ! Dirichelet boundary condition at the east of the domain
2915  DO i=1,n_vars
2916 
2917  IF ( bce(i)%flag .EQ. 0 ) THEN
2918 
2919  qrece(i) = bce(i)%value
2920 
2921  END IF
2922 
2923  ENDDO
2924 
2925  DO i=n_vars+1,n_vars+2
2926 
2927  CALL qp_to_qp2( qrecw(1:n_vars) , b_cent(j,k) , qp2recw )
2928  qrecw(i) = qp2recw(i-n_vars+1)
2929 
2930  CALL qp_to_qp2( qrece(1:n_vars) , b_cent(j,k) , qp2rece )
2931  qrece(i) = qp2rece(i-n_vars+1)
2932 
2933  END DO
2934 
2935  ELSE
2936 
2937  ! correction for radial source inlet x-interfaces:
2938  ! the physical variables at the x-interfaces qrecW or
2939  ! qrecE are computed from the radial inlet values
2940  IF ( radial_source_flag .AND. ( source_cell(j,k).EQ.2 ) ) THEN
2941 
2942  IF ( sourcee(j,k) ) THEN
2943 
2944  CALL eval_source_bdry( t, sourcee_vect_x(j,k) , &
2945  sourcee_vect_y(j,k) , source_bdry )
2946 
2947  qrece(1:n_vars+2) = source_bdry(1:n_vars+2)
2948 
2949  ELSEIF ( sourcew(j,k) ) THEN
2950 
2951  CALL eval_source_bdry( t, sourcew_vect_x(j,k) , &
2952  sourcew_vect_y(j,k) , source_bdry )
2953 
2954  qrecw(1:n_vars+2) = source_bdry(1:n_vars+2)
2955 
2956  END IF
2957 
2958  END IF
2959 
2960  END IF
2961 
2962  CALL qp_to_qc( qrecw,b_interfacer(j,k),q_interfacer(:,j,k) )
2963  CALL qp_to_qc( qrece,b_interfacel(j+1,k),q_interfacel(:,j+1,k) )
2964 
2965  qp_interfacer(1:n_vars+2,j,k) = qrecw(1:n_vars+2)
2966  qp_interfacel(1:n_vars+2,j+1,k) = qrece(1:n_vars+2)
2967 
2968  IF ( j.EQ.1 ) THEN
2969 
2970  ! Interface value at the left of first x-interface (external)
2971  q_interfacel(:,j,k) = q_interfacer(:,j,k)
2972  qp_interfacel(:,j,k) = qp_interfacer(:,j,k)
2973 
2974  !WRITE(*,*) 'q_interfaceL(:,j,k)',q_interfaceL(:,j,k)
2975  !READ(*,*)
2976 
2977  ELSEIF ( j.EQ.comp_cells_x ) THEN
2978 
2979  ! Interface value at the right of last x-interface (external)
2980  q_interfacer(:,j+1,k) = q_interfacel(:,j+1,k)
2981  qp_interfacer(:,j+1,k) = qp_interfacel(:,j+1,k)
2982 
2983  ELSE
2984 
2985  IF ( radial_source_flag .AND. ( source_cell(j,k) .EQ. 2 ) ) THEN
2986 
2987  IF ( sourcee(j,k) ) THEN
2988 
2989  q_interfacer(:,j+1,k) = q_interfacel(:,j+1,k)
2990  qp_interfacer(:,j+1,k) = qp_interfacel(:,j+1,k)
2991 
2992  ELSEIF ( sourcew(j,k) ) THEN
2993 
2994  q_interfacel(:,j,k) = q_interfacer(:,j,k)
2995  qp_interfacel(:,j,k) = qp_interfacer(:,j,k)
2996 
2997  END IF
2998 
2999  END IF
3000 
3001  END IF
3002 
3003  ELSE
3004 
3005  ! for case comp_cells_x = 1
3006  q_interfacer(1:n_vars,j,k) = q_expl(1:n_vars,j,k)
3007  q_interfacel(1:n_vars,j+1,k) = q_expl(1:n_vars,j,k)
3008 
3009  qp_interfacer(1:n_vars+2,j,k) = qrec(1:n_vars+2,j,k)
3010  qp_interfacel(1:n_vars+2,j+1,k) = qrec(1:n_vars+2,j,k)
3011 
3012  END IF
3013 
3014  IF ( comp_cells_y .EQ. 1 ) THEN
3015 
3016  ! comp_cells_y = 1
3017  q_interfacet(:,j,k) = q_expl(:,j,k)
3018  q_interfaceb(:,j,k+1) = q_expl(:,j,k)
3019 
3020  qp_interfacet(:,j,k) = qrec(:,j,k)
3021  qp_interfaceb(:,j,k+1) = qrec(:,j,k)
3022 
3023  ELSE
3024 
3025  IF ( ( k .GT. 1 ) .AND. ( k .LT. comp_cells_y ) ) THEN
3026 
3027  IF ( q_expl(1,j,k) .EQ. 0.d0 ) THEN
3028 
3029  IF ( ( .NOT. radial_source_flag ) .OR. &
3030  ( ( radial_source_flag ) .AND. &
3031  ( source_cell(j,k) .EQ. 0 ) ) ) THEN
3032 
3033  ! In the internal cell, if thickness h is 0 at the center
3034  ! of the cell, then all the variables are 0 at the center
3035  ! and at the interfaces (no conversion back is needed from
3036  ! reconstructed to conservative)
3037 
3038  q_interfacet(:,j,k) = 0.d0
3039  q_interfaceb(:,j,k+1) = 0.d0
3040 
3041  qp_interfacet(1:3,j,k) = 0.d0
3042  qp_interfacet(4:n_vars,j,k) = qrecs(4:n_vars)
3043  qp_interfacet(n_vars+1:n_vars+2,j,k) = 0.d0
3044 
3045  qp_interfaceb(1:3,j,k+1) = 0.d0
3046  qp_interfaceb(4:n_vars,j,k+1) = qrecn(4:n_vars)
3047  qp_interfaceb(n_vars+1:n_vars+2,j,k+1) = 0.d0
3048 
3049  END IF
3050 
3051  END IF
3052 
3053  END IF
3054 
3055  IF ( k .EQ. 1 ) THEN
3056 
3057  ! Dirichelet boundary condition at the south of the domain
3058  DO i=1,n_vars
3059 
3060  IF ( bcs(i)%flag .EQ. 0 ) THEN
3061 
3062  qrecs(i) = bcs(i)%value
3063 
3064  END IF
3065 
3066  ENDDO
3067 
3068  DO i=n_vars+1,n_vars+2
3069 
3070  CALL qp_to_qp2( qrecs(1:n_vars) , b_cent(j,k) , qp2recs )
3071  qrecs(i) = qp2recs(i-n_vars+1)
3072 
3073  CALL qp_to_qp2( qrecn(1:n_vars) , b_cent(j,k) , qp2recn )
3074  qrecn(i) = qp2recn(i-n_vars+1)
3075 
3076  END DO
3077 
3078  ELSEIF ( k .EQ. comp_cells_y ) THEN
3079 
3080  ! Dirichelet boundary condition at the north of the domain
3081  DO i=1,n_vars
3082 
3083  IF ( bcn(i)%flag .EQ. 0 ) THEN
3084 
3085  qrecn(i) = bcn(i)%value
3086 
3087  END IF
3088 
3089  ENDDO
3090 
3091  DO i=n_vars+1,n_vars+2
3092 
3093  CALL qp_to_qp2( qrecs(1:n_vars) , b_cent(j,k) , qp2recs )
3094  qrecs(i) = qp2recs(i-n_vars+1)
3095 
3096  CALL qp_to_qp2( qrecn(1:n_vars) , b_cent(j,k) , qp2recn )
3097  qrecn(i) = qp2recn(i-n_vars+1)
3098 
3099  END DO
3100 
3101  ELSE
3102 
3103  IF ( radial_source_flag .AND. ( source_cell(j,k) .EQ. 2 ) ) THEN
3104 
3105  IF ( sources(j,k) ) THEN
3106 
3107  CALL eval_source_bdry( t, sources_vect_x(j,k) , &
3108  sources_vect_y(j,k) , source_bdry )
3109 
3110  qrecs(1:n_vars+2) = source_bdry(1:n_vars+2)
3111 
3112  ELSEIF ( sourcen(j,k) ) THEN
3113 
3114  CALL eval_source_bdry( t, sourcen_vect_x(j,k) , &
3115  sourcen_vect_y(j,k) , source_bdry )
3116 
3117  qrecn(1:n_vars+2) = source_bdry(1:n_vars+2)
3118 
3119  END IF
3120 
3121  END IF
3122 
3123  END IF
3124 
3125  CALL qp_to_qc( qrecs, b_interfacet(j,k), q_interfacet(:,j,k))
3126  CALL qp_to_qc( qrecn, b_interfaceb(j,k+1), q_interfaceb(:,j,k+1))
3127 
3128  qp_interfacet(1:n_vars+2,j,k) = qrecs(1:n_vars+2)
3129  qp_interfaceb(1:n_vars+2,j,k+1) = qrecn(1:n_vars+2)
3130 
3131  IF ( k .EQ. 1 ) THEN
3132 
3133  ! Interface value at the bottom of first y-interface (external)
3134  q_interfaceb(:,j,k) = q_interfacet(:,j,k)
3135  qp_interfaceb(:,j,k) = qp_interfacet(:,j,k)
3136 
3137  ELSEIF ( k .EQ. comp_cells_y ) THEN
3138 
3139  ! Interface value at the top of last y-interface (external)
3140  q_interfacet(:,j,k+1) = q_interfaceb(:,j,k+1)
3141  qp_interfacet(:,j,k+1) = qp_interfaceb(:,j,k+1)
3142 
3143  ELSE
3144 
3145  IF ( radial_source_flag .AND. ( source_cell(j,k) .EQ. 2 ) ) THEN
3146 
3147  IF ( sources(j,k) ) THEN
3148 
3149  q_interfaceb(:,j,k) = q_interfacet(:,j,k)
3150  qp_interfaceb(:,j,k) = qp_interfacet(:,j,k)
3151 
3152  ELSEIF ( sourcen(j,k) ) THEN
3153 
3154  q_interfacet(:,j,k+1) = q_interfaceb(:,j,k+1)
3155  qp_interfacet(:,j,k+1) = qp_interfaceb(:,j,k+1)
3156 
3157  END IF
3158 
3159  END IF
3160 
3161  END IF
3162 
3163  END IF
3164 
3165  END DO
3166 
3167  DEALLOCATE ( qrec )
3168  DEALLOCATE ( qrecw )
3169  DEALLOCATE ( qrece )
3170  DEALLOCATE ( qrecs )
3171  DEALLOCATE ( qrecn )
3172  DEALLOCATE ( source_bdry )
3173  DEALLOCATE ( qrec_prime_x )
3174  DEALLOCATE ( qrec_prime_y )
3175 
3176  RETURN
3177 
3178  END SUBROUTINE reconstruction
3179 
3180 
3181  !******************************************************************************
3183  !
3189  !******************************************************************************
3190 
3191  SUBROUTINE eval_speeds
3192 
3193  ! External procedures
3194  USE constitutive_2d, ONLY : eval_local_speeds_x, eval_local_speeds_y
3195 
3196  IMPLICIT NONE
3197 
3198  REAL*8 :: abslambdaL_min(n_vars) , abslambdaL_max(n_vars)
3199  REAL*8 :: abslambdaR_min(n_vars) , abslambdaR_max(n_vars)
3200  REAL*8 :: abslambdaB_min(n_vars) , abslambdaB_max(n_vars)
3201  REAL*8 :: abslambdaT_min(n_vars) , abslambdaT_max(n_vars)
3202  REAL*8 :: min_r(n_vars) , max_r(n_vars)
3203 
3204  INTEGER :: j,k,l
3205 
3206  IF ( comp_cells_x .GT. 1 ) THEN
3207 
3208  x_interfaces_loop:DO l = 1,solve_interfaces_x
3209 
3210  j = j_stag_x(l)
3211  k = k_stag_x(l)
3212 
3213  CALL eval_local_speeds_x( qp_interfacel(:,j,k) , abslambdal_min , &
3214  abslambdal_max )
3215 
3216  CALL eval_local_speeds_x( qp_interfacer(:,j,k) , abslambdar_min , &
3217  abslambdar_max )
3218 
3219  min_r = min(abslambdal_min , abslambdar_min , 0.0d0)
3220  max_r = max(abslambdal_max , abslambdar_max , 0.0d0)
3221 
3222  a_interface_xneg(:,j,k) = min_r
3223  a_interface_xpos(:,j,k) = max_r
3224 
3225  END DO x_interfaces_loop
3226 
3227  END IF
3228 
3229  IF ( comp_cells_y .GT. 1 ) THEN
3230 
3231  y_interfaces_loop:DO l = 1,solve_interfaces_y
3232 
3233  j = j_stag_y(l)
3234  k = k_stag_y(l)
3235 
3236  CALL eval_local_speeds_y( qp_interfaceb(:,j,k) , abslambdab_min , &
3237  abslambdab_max )
3238 
3239  CALL eval_local_speeds_y( qp_interfacet(:,j,k) , abslambdat_min , &
3240  abslambdat_max )
3241 
3242  min_r = min(abslambdab_min , abslambdat_min , 0.0d0)
3243  max_r = max(abslambdab_max , abslambdat_max , 0.0d0)
3244 
3245  a_interface_yneg(:,j,k) = min_r
3246  a_interface_ypos(:,j,k) = max_r
3247 
3248  END DO y_interfaces_loop
3249 
3250  END IF
3251 
3252  RETURN
3253 
3254  END SUBROUTINE eval_speeds
3255 
3256 END MODULE solver_2d
geometry_2d::b_prime_x
real *8, dimension(:,:), allocatable b_prime_x
Topography slope (x direction) at the centers of the control volumes.
Definition: geometry_2d.f90:44
solver_2d::q_cellne
real *8, dimension(:,:,:), allocatable q_cellne
Reconstructed value at the NE corner of cell.
Definition: solver_2d.f90:69
parameters_2d::max_dt
real *8 max_dt
Largest time step allowed.
Definition: parameters_2d.f90:16
solver_2d::eval_speeds
subroutine eval_speeds
Characteristic speeds.
Definition: solver_2d.f90:3192
geometry_2d::b_interfacel
real *8, dimension(:,:), allocatable b_interfacel
Reconstructed value at the left of the x-interface.
Definition: geometry_2d.f90:29
solver_2d::a_dirk
real *8, dimension(:), allocatable a_dirk
Explicit coeff. for the non-hyp. part for a single step of the R-K scheme.
Definition: solver_2d.f90:142
geometry_2d::dy
real *8 dy
Control volumes size.
Definition: geometry_2d.f90:93
solver_2d::qp_cellnw
real *8, dimension(:,:,:), allocatable qp_cellnw
Reconstructed physical value at the NW corner of cell.
Definition: solver_2d.f90:76
solver_2d::nhj
real *8, dimension(:,:), allocatable nhj
Local Intermediate non-hyperbolic terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:170
solver_2d::eval_explicit_terms
subroutine eval_explicit_terms(qc_expl, qp_expl, expl_terms)
Evaluate the explicit terms.
Definition: solver_2d.f90:1999
parameters_2d::n_rk
integer n_rk
Runge-Kutta order.
Definition: parameters_2d.f90:131
parameters_2d::tol_rel
real *8, parameter tol_rel
Definition: parameters_2d.f90:136
solver_2d::q_rk
real *8, dimension(:,:,:,:), allocatable q_rk
Intermediate solutions of the Runge-Kutta scheme.
Definition: solver_2d.f90:148
solver_2d::j_stag_y
integer, dimension(:), allocatable j_stag_y
Definition: solver_2d.f90:196
geometry_2d::x_comp
real *8, dimension(:), allocatable x_comp
Location of the centers (x) of the control volume of the domain.
Definition: geometry_2d.f90:14
solver_2d::q0
real *8, dimension(:,:,:), allocatable q0
Conservative variables.
Definition: solver_2d.f90:44
geometry_2d::b_cent
real *8, dimension(:,:), allocatable b_cent
Topography at the centers of the control volumes.
Definition: geometry_2d.f90:41
solver_2d::reconstruction
subroutine reconstruction(q_expl, qp_expl)
Linear reconstruction.
Definition: solver_2d.f90:2486
geometry_2d::comp_cells_x
integer comp_cells_x
Number of control volumes x in the comp. domain.
Definition: geometry_2d.f90:98
solver_2d::eval_flux_gforce
subroutine eval_flux_gforce
Numerical fluxes GFORCE.
Definition: solver_2d.f90:2449
constitutive_2d::eval_local_speeds_x
subroutine eval_local_speeds_x(qpj, vel_min, vel_max)
Local Characteristic speeds x direction.
Definition: constitutive_2d.f90:973
geometry_2d::sourcen_vect_y
real *8, dimension(:,:), allocatable sourcen_vect_y
Definition: geometry_2d.f90:84
geometry_2d::grav_surf
real *8, dimension(:,:), allocatable grav_surf
gravity vector wrt surface coordinates for each cell
Definition: geometry_2d.f90:53
solver_2d::divflux
real *8, dimension(:,:,:,:), allocatable divflux
Intermediate hyperbolic terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:155
geometry_2d::sourcee_vect_x
real *8, dimension(:,:), allocatable sourcee_vect_x
Definition: geometry_2d.f90:74
geometry_2d::y_comp
real *8, dimension(:), allocatable y_comp
Location of the centers (y) of the control volume of the domain.
Definition: geometry_2d.f90:20
solver_2d::mask11
logical, dimension(:,:), allocatable mask11
Definition: solver_2d.f90:123
solver_2d::source_xy
real *8, dimension(:,:), allocatable source_xy
Array defining fraction of cells affected by source term.
Definition: solver_2d.f90:110
constitutive_2d::mixt_var
subroutine mixt_var(qpj, r_Ri, r_rho_m, r_rho_c, r_red_grav)
Physical variables.
Definition: constitutive_2d.f90:569
solver_2d::k_cent
integer, dimension(:), allocatable k_cent
Definition: solver_2d.f90:191
solver_2d::average_kt
subroutine average_kt(a1, a2, w1, w2, w_avg)
averaged KT flux
Definition: solver_2d.f90:2410
parameters_2d
Parameters.
Definition: parameters_2d.f90:7
geometry_2d::dx
real *8 dx
Control volumes size.
Definition: geometry_2d.f90:90
solver_2d::si_nh
real *8, dimension(:,:,:,:), allocatable si_nh
Intermediate semi-implicit non-hyperbolic terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:161
solver_2d::q_cellse
real *8, dimension(:,:,:), allocatable q_cellse
Reconstructed value at the SE corner of cell.
Definition: solver_2d.f90:73
solver_2d::qp_interfaceb
real *8, dimension(:,:,:), allocatable qp_interfaceb
Reconstructed physical value at the bottom of the y-interface.
Definition: solver_2d.f90:62
geometry_2d::comp_cells_y
integer comp_cells_y
Number of control volumes y in the comp. domain.
Definition: geometry_2d.f90:100
solver_2d
Numerical solver.
Definition: solver_2d.f90:12
solver_2d::a_interface_xneg
real *8, dimension(:,:,:), allocatable a_interface_xneg
Local speeds at the left of the x-interface.
Definition: solver_2d.f90:95
solver_2d::deallocate_solver_variables
subroutine deallocate_solver_variables
Memory deallocation.
Definition: solver_2d.f90:452
constitutive_2d::qp_to_qp2
subroutine qp_to_qp2(qpj, Bj, qp2j)
Additional Physical variables.
Definition: constitutive_2d.f90:934
geometry_2d::sourcee
logical, dimension(:,:), allocatable sourcee
Definition: geometry_2d.f90:64
solver_2d::expl_terms
real *8, dimension(:,:,:,:), allocatable expl_terms
Intermediate explicit terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:164
geometry_2d::sourcew
logical, dimension(:,:), allocatable sourcew
Definition: geometry_2d.f90:65
parameters_2d::bcw
type(bc), dimension(:), allocatable bcw
bcW&flag defines the west boundary condition:
Definition: parameters_2d.f90:178
solver_2d::nh
real *8, dimension(:,:,:,:), allocatable nh
Intermediate non-hyperbolic terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:158
geometry_2d::sourcew_vect_y
real *8, dimension(:,:), allocatable sourcew_vect_y
Definition: geometry_2d.f90:78
geometry_2d::b_interfacer
real *8, dimension(:,:), allocatable b_interfacer
Reconstructed value at the right of the x-interface.
Definition: geometry_2d.f90:32
geometry_2d::minmod
real *8 function minmod(a, b)
Definition: geometry_2d.f90:1014
constitutive_2d::eval_nh_semi_impl_terms
subroutine eval_nh_semi_impl_terms(grav3_surf, qcj, nh_semi_impl_term)
Non-Hyperbolic semi-implicit terms.
Definition: constitutive_2d.f90:1478
solver_2d::k_stag_x
integer, dimension(:), allocatable k_stag_x
Definition: solver_2d.f90:194
solver_2d::qp_interfacet
real *8, dimension(:,:,:), allocatable qp_interfacet
Reconstructed physical value at the top of the y-interface.
Definition: solver_2d.f90:64
geometry_2d::sourcew_vect_x
real *8, dimension(:,:), allocatable sourcew_vect_x
Definition: geometry_2d.f90:77
solver_2d::q_cellsw
real *8, dimension(:,:,:), allocatable q_cellsw
Reconstructed value at the SW corner of cell.
Definition: solver_2d.f90:71
solver_2d::a_tilde
real *8, dimension(:), allocatable a_tilde
Explicit coeff. for the hyperbolic part for a single step of the R-K scheme.
Definition: solver_2d.f90:139
solver_2d::dt
real *8 dt
Time step.
Definition: solver_2d.f90:121
parameters_2d::n_vars
integer n_vars
Number of conservative variables.
Definition: parameters_2d.f90:124
parameters_2d::topo_change_flag
logical topo_change_flag
Definition: parameters_2d.f90:71
geometry_2d::b_interfaceb
real *8, dimension(:,:), allocatable b_interfaceb
Reconstructed value at the bottom of the y-interface.
Definition: geometry_2d.f90:35
solver_2d::a_interface_yneg
real *8, dimension(:,:,:), allocatable a_interface_yneg
Local speeds at the bottom of the y-interface.
Definition: solver_2d.f90:99
solver_2d::eval_flux_lxf
subroutine eval_flux_lxf
Numerical fluxes Lax-Friedrichs.
Definition: solver_2d.f90:2463
solver_2d::qp_cellne
real *8, dimension(:,:,:), allocatable qp_cellne
Reconstructed physical value at the NE corner of cell.
Definition: solver_2d.f90:78
geometry_2d::comp_interfaces_y
integer comp_interfaces_y
Number of interfaces (comp_cells_y+1)
Definition: geometry_2d.f90:101
parameters_2d::cfl
real *8 cfl
Courant-Friedrichs-Lewy parameter.
Definition: parameters_2d.f90:18
constitutive_2d::eval_expl_terms
subroutine eval_expl_terms(Bprimej_x, Bprimej_y, source_xy, qpj, qcj, expl_term)
Explicit Forces term.
Definition: constitutive_2d.f90:1604
parameters_2d::max_nl_iter
integer, parameter max_nl_iter
Definition: parameters_2d.f90:133
solver_2d::residual_term
real *8, dimension(:,:,:), allocatable residual_term
Sum of all the terms of the equations except the transient term.
Definition: solver_2d.f90:188
constitutive_2d::eval_fluxes
subroutine eval_fluxes(qcj, qpj, dir, flux)
Hyperbolic Fluxes.
Definition: constitutive_2d.f90:1077
solver_2d::update_erosion_deposition_cell
subroutine update_erosion_deposition_cell(dt)
Evaluate the eroion/deposition terms.
Definition: solver_2d.f90:1880
geometry_2d::comp_cells_xy
integer comp_cells_xy
Definition: geometry_2d.f90:103
constitutive_2d
Constitutive equations.
Definition: constitutive_2d.f90:4
solver_2d::normalize_q
logical normalize_q
Flag for the normalization of the array q in the implicit solution scheme.
Definition: solver_2d.f90:179
solver_2d::mask21
logical, dimension(:,:), allocatable mask21
Definition: solver_2d.f90:123
solver_2d::h_interface_x
real *8, dimension(:,:,:), allocatable h_interface_x
Semidiscrete numerical interface fluxes.
Definition: solver_2d.f90:103
solver_2d::a_dirk_ij
real *8, dimension(:,:), allocatable a_dirk_ij
Butcher Tableau for the implicit part of the Runge-Kutta scheme.
Definition: solver_2d.f90:130
solver_2d::solve_mask_y
logical, dimension(:,:), allocatable solve_mask_y
Definition: solver_2d.f90:114
solver_2d::a_interface_ypos
real *8, dimension(:,:,:), allocatable a_interface_ypos
Local speeds at the top of the y-interface.
Definition: solver_2d.f90:101
constitutive_2d::implicit_flag
logical, dimension(:), allocatable implicit_flag
flag used for size of implicit non linear-system
Definition: constitutive_2d.f90:14
solver_2d::h_interface_y
real *8, dimension(:,:,:), allocatable h_interface_y
Semidiscrete numerical interface fluxes.
Definition: solver_2d.f90:105
geometry_2d::limit
subroutine limit(v, z, limiter, slope_lim)
Slope limiter.
Definition: geometry_2d.f90:946
solver_2d::qp_interfacer
real *8, dimension(:,:,:), allocatable qp_interfacer
Reconstructed physical value at the right of the x-interface.
Definition: solver_2d.f90:60
constitutive_2d::erosion_coeff
real *8, dimension(:), allocatable erosion_coeff
erosion model coefficient (units: m-1 )
Definition: constitutive_2d.f90:139
constitutive_2d::settling_flag
logical settling_flag
Flag to determine if sedimentation is active.
Definition: constitutive_2d.f90:133
solver_2d::q_interfacel
real *8, dimension(:,:,:), allocatable q_interfacel
Reconstructed value at the left of the x-interface.
Definition: solver_2d.f90:49
solver_2d::solve_rk_step
subroutine solve_rk_step(qj, qj_old, a_tilde, a_dirk, a_diag)
Runge-Kutta single step integration.
Definition: solver_2d.f90:1228
constitutive_2d::eval_source_bdry
subroutine eval_source_bdry(time, vect_x, vect_y, source_bdry)
Internal boundary source fluxes.
Definition: constitutive_2d.f90:1888
geometry_2d::sourcee_vect_y
real *8, dimension(:,:), allocatable sourcee_vect_y
Definition: geometry_2d.f90:75
geometry_2d::sources_vect_y
real *8, dimension(:,:), allocatable sources_vect_y
Definition: geometry_2d.f90:81
solver_2d::omega_tilde
real *8, dimension(:), allocatable omega_tilde
Coefficients for the explicit part of the Runge-Kutta scheme.
Definition: solver_2d.f90:133
constitutive_2d::eval_local_speeds_y
subroutine eval_local_speeds_y(qpj, vel_min, vel_max)
Local Characteristic speeds y direction.
Definition: constitutive_2d.f90:1024
parameters_2d::solver_scheme
character(len=20) solver_scheme
Finite volume method: .
Definition: parameters_2d.f90:151
geometry_2d::sourcen_vect_x
real *8, dimension(:,:), allocatable sourcen_vect_x
Definition: geometry_2d.f90:83
geometry_2d::source_cell
integer, dimension(:,:), allocatable source_cell
Definition: geometry_2d.f90:63
constitutive_2d::eval_erosion_dep_term
subroutine eval_erosion_dep_term(qpj, dt, erosion_term, deposition_term)
Erosion/Deposition term.
Definition: constitutive_2d.f90:1668
geometry_2d
Grid module.
Definition: geometry_2d.f90:7
solver_2d::omega
real *8, dimension(:), allocatable omega
Coefficients for the implicit part of the Runge-Kutta scheme.
Definition: solver_2d.f90:136
solver_2d::expl_terms_j
real *8, dimension(:,:), allocatable expl_terms_j
Local Intermediate explicit terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:173
solver_2d::t
real *8 t
time
Definition: solver_2d.f90:39
solver_2d::q_cellnw
real *8, dimension(:,:,:), allocatable q_cellnw
Reconstructed value at the NW corner of cell.
Definition: solver_2d.f90:67
solver_2d::check_solve
subroutine check_solve
Masking of cells to solve.
Definition: solver_2d.f90:548
solver_2d::qp_cellsw
real *8, dimension(:,:,:), allocatable qp_cellsw
Reconstructed physical value at the SW corner of cell.
Definition: solver_2d.f90:80
solver_2d::qp_cellse
real *8, dimension(:,:,:), allocatable qp_cellse
Reconstructed physical value at the SE corner of cell.
Definition: solver_2d.f90:82
solver_2d::a_interface_xpos
real *8, dimension(:,:,:), allocatable a_interface_xpos
Local speeds at the right of the x-interface.
Definition: solver_2d.f90:97
solver_2d::k_stag_y
integer, dimension(:), allocatable k_stag_y
Definition: solver_2d.f90:197
solver_2d::mask22
logical, dimension(:,:), allocatable mask22
Definition: solver_2d.f90:123
solver_2d::lnsrch
subroutine lnsrch(qj_rel_NR_old, qj_org, qj_old, scal_f_old, grad_f, desc_dir, coeff_f, qj_rel, scal_f, right_term, stpmax, check)
Search the descent stepsize.
Definition: solver_2d.f90:1543
parameters_2d::verbose_level
integer verbose_level
Definition: parameters_2d.f90:161
solver_2d::eval_hyperbolic_terms
subroutine eval_hyperbolic_terms(q_expl, qp_expl, divFlux)
Semidiscrete finite volume central scheme.
Definition: solver_2d.f90:2053
geometry_2d::b_prime_y
real *8, dimension(:,:), allocatable b_prime_y
Topography slope (y direction) at the centers of the control volumes.
Definition: geometry_2d.f90:47
solver_2d::timestep
subroutine timestep
Time-step computation.
Definition: solver_2d.f90:697
solver_2d::si_nhj
real *8, dimension(:,:), allocatable si_nhj
Local Intermediate semi-impl non-hyperbolic terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:176
solver_2d::a_interface_x_max
real *8, dimension(:,:,:), allocatable a_interface_x_max
Max local speeds at the x-interface.
Definition: solver_2d.f90:89
solver_2d::solve_interfaces_x
integer solve_interfaces_x
Definition: solver_2d.f90:117
constitutive_2d::qc_to_qp
subroutine qc_to_qp(qc, qp)
Conservative to physical variables.
Definition: constitutive_2d.f90:678
solver_2d::q_interfaceb
real *8, dimension(:,:,:), allocatable q_interfaceb
Reconstructed value at the bottom of the y-interface.
Definition: solver_2d.f90:53
solver_2d::j_cent
integer, dimension(:), allocatable j_cent
Definition: solver_2d.f90:190
geometry_2d::x_stag
real *8, dimension(:), allocatable x_stag
Location of the boundaries (x) of the control volumes of the domain.
Definition: geometry_2d.f90:17
geometry_2d::sources_vect_x
real *8, dimension(:,:), allocatable sources_vect_x
Definition: geometry_2d.f90:80
solver_2d::a_tilde_ij
real *8, dimension(:,:), allocatable a_tilde_ij
Butcher Tableau for the explicit part of the Runge-Kutta scheme.
Definition: solver_2d.f90:128
solver_2d::solve_interfaces_y
integer solve_interfaces_y
Definition: solver_2d.f90:118
solver_2d::a_diag
real *8 a_diag
Implicit coeff. for the non-hyp. part for a single step of the R-K scheme.
Definition: solver_2d.f90:145
solver_2d::opt_search_nl
logical opt_search_nl
Flag for the search of optimal step size in the implicit solution scheme.
Definition: solver_2d.f90:185
solver_2d::q_fv
real *8, dimension(:,:,:), allocatable q_fv
Solution of the finite-volume semidiscrete cheme.
Definition: solver_2d.f90:46
solver_2d::normalize_f
logical normalize_f
Flag for the normalization of the array f in the implicit solution scheme.
Definition: solver_2d.f90:182
parameters_2d::limiter
integer, dimension(10) limiter
Limiter for the slope in the linear reconstruction: .
Definition: parameters_2d.f90:144
solver_2d::qp_rk
real *8, dimension(:,:,:,:), allocatable qp_rk
Intermediate physical solutions of the Runge-Kutta scheme.
Definition: solver_2d.f90:151
parameters_2d::bcn
type(bc), dimension(:), allocatable bcn
bcN&flag defines the north boundary condition:
Definition: parameters_2d.f90:208
geometry_2d::b_interfacet
real *8, dimension(:,:), allocatable b_interfacet
Reconstructed value at the top of the y-interface.
Definition: geometry_2d.f90:38
geometry_2d::sourcen
logical, dimension(:,:), allocatable sourcen
Definition: geometry_2d.f90:67
geometry_2d::dx2
real *8 dx2
Half x Control volumes size.
Definition: geometry_2d.f90:96
parameters_2d::n_eqns
integer n_eqns
Number of equations.
Definition: parameters_2d.f90:125
solver_2d::eval_jacobian
subroutine eval_jacobian(qj_rel, qj_org, coeff_f, left_matrix)
Evaluate the jacobian.
Definition: solver_2d.f90:1810
constitutive_2d::qp_to_qc
subroutine qp_to_qc(qp, B, qc)
Physical to conservative variables.
Definition: constitutive_2d.f90:752
geometry_2d::deposit
real *8, dimension(:,:,:), allocatable deposit
deposit for the different classes
Definition: geometry_2d.f90:59
parameters_2d::tol_abs
real *8, parameter tol_abs
Definition: parameters_2d.f90:135
constitutive_2d::eval_nonhyperbolic_terms
subroutine eval_nonhyperbolic_terms(c_qj, c_nh_term_impl, r_qj, r_nh_term_impl)
Non-Hyperbolic terms.
Definition: constitutive_2d.f90:1208
solver_2d::eval_flux_up
subroutine eval_flux_up
Upwind numerical fluxes.
Definition: solver_2d.f90:2157
solver_2d::t_imex1
real *8 t_imex1
Definition: solver_2d.f90:199
parameters_2d::reconstr_coeff
real *8 reconstr_coeff
Slope coefficient in the linear reconstruction.
Definition: parameters_2d.f90:22
solver_2d::imex_rk_solver
subroutine imex_rk_solver
Runge-Kutta integration.
Definition: solver_2d.f90:828
parameters_2d::bcs
type(bc), dimension(:), allocatable bcs
bcS&flag defines the south boundary condition:
Definition: parameters_2d.f90:198
solver_2d::solve_mask
logical, dimension(:,:), allocatable solve_mask
Definition: solver_2d.f90:112
geometry_2d::comp_interfaces_x
integer comp_interfaces_x
Number of interfaces (comp_cells_x+1)
Definition: geometry_2d.f90:99
geometry_2d::y_stag
real *8, dimension(:), allocatable y_stag
Location of the boundaries (x) of the control volumes of the domain.
Definition: geometry_2d.f90:23
parameters_2d::bce
type(bc), dimension(:), allocatable bce
bcE&flag defines the east boundary condition:
Definition: parameters_2d.f90:188
solver_2d::eval_f
subroutine eval_f(qj, qj_old, a_tilde, a_dirk, a_diag, coeff_f, f_nl, scal_f)
Evaluate the nonlinear system.
Definition: solver_2d.f90:1759
solver_2d::solve_mask_x
logical, dimension(:,:), allocatable solve_mask_x
Definition: solver_2d.f90:113
solver_2d::q1max
real *8, dimension(:,:), allocatable q1max
Maximum over time of thickness.
Definition: solver_2d.f90:86
solver_2d::q_interfacer
real *8, dimension(:,:,:), allocatable q_interfacer
Reconstructed value at the right of the x-interface.
Definition: solver_2d.f90:51
parameters_2d::radial_source_flag
logical radial_source_flag
Definition: parameters_2d.f90:73
solver_2d::j_stag_x
integer, dimension(:), allocatable j_stag_x
Definition: solver_2d.f90:193
solver_2d::divfluxj
real *8, dimension(:,:), allocatable divfluxj
Local Intermediate hyperbolic terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:167
solver_2d::eval_flux_kt
subroutine eval_flux_kt
Semidiscrete numerical fluxes.
Definition: solver_2d.f90:2275
solver_2d::solve_cells
integer solve_cells
Definition: solver_2d.f90:116
solver_2d::i_rk
integer i_rk
loop counter for the RK iteration
Definition: solver_2d.f90:125
solver_2d::allocate_solver_variables
subroutine allocate_solver_variables
Memory allocation.
Definition: solver_2d.f90:216
solver_2d::t_imex2
real *8 t_imex2
Definition: solver_2d.f90:199
geometry_2d::dy2
real *8 dy2
Half y Control volumes size.
Definition: geometry_2d.f90:97
solver_2d::mask12
logical, dimension(:,:), allocatable mask12
Definition: solver_2d.f90:123
solver_2d::a_interface_y_max
real *8, dimension(:,:,:), allocatable a_interface_y_max
Max local speeds at the y-interface.
Definition: solver_2d.f90:91
solver_2d::qp_interfacel
real *8, dimension(:,:,:), allocatable qp_interfacel
Reconstructed physical value at the left of the x-interface.
Definition: solver_2d.f90:58
solver_2d::q_interfacet
real *8, dimension(:,:,:), allocatable q_interfacet
Reconstructed value at the top of the y-interface.
Definition: solver_2d.f90:55
parameters_2d::n_solid
integer n_solid
Number of solid classes.
Definition: parameters_2d.f90:127
constitutive_2d::rho_s
real *8, dimension(:), allocatable rho_s
Density of sediments ( units: kg m-3 )
Definition: constitutive_2d.f90:124
solver_2d::qp
real *8, dimension(:,:,:), allocatable qp
Physical variables ( )
Definition: solver_2d.f90:107
constitutive_2d::eval_topo_term
subroutine eval_topo_term(qpj, deposition_avg_term, erosion_avg_term, eqns_term, deposit_term)
Topography modification related term.
Definition: constitutive_2d.f90:1774
parameters_2d::n_nh
integer n_nh
Number of non-hyperbolic terms.
Definition: parameters_2d.f90:129
geometry_2d::sources
logical, dimension(:,:), allocatable sources
Definition: geometry_2d.f90:66
solver_2d::q
real *8, dimension(:,:,:), allocatable q
Conservative variables.
Definition: solver_2d.f90:42
Generated on Fri Dec 13 2019 21:38:13 for SW_VAR_DENS_MODEL by   doxygen 1.8.11