LAHARS-MODEL  0.1
templategithubproject
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
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 , b_stag_x , b_stag_y
22  USE geometry_2d, ONLY : b_ver
23  USE geometry_2d, ONLY : grav_surf
24 
25  USE parameters_2d, ONLY : n_eqns , n_vars , n_nh
26  USE parameters_2d, ONLY : n_rk
27  USE parameters_2d, ONLY : verbose_level
28  USE parameters_2d, ONLY : solid_transport_flag
29 
30  USE parameters_2d, ONLY : bcw , bce , bcs , bcn
31 
32  IMPLICIT none
33 
35  REAL*8, ALLOCATABLE :: q(:,:,:)
37  REAL*8, ALLOCATABLE :: q0(:,:,:)
39  REAL*8, ALLOCATABLE :: q_fv(:,:,:)
40 
42  REAL*8, ALLOCATABLE :: q_interfacel(:,:,:)
44  REAL*8, ALLOCATABLE :: q_interfacer(:,:,:)
46  REAL*8, ALLOCATABLE :: q_interfaceb(:,:,:)
48  REAL*8, ALLOCATABLE :: q_interfacet(:,:,:)
49 
51  REAL*8, ALLOCATABLE :: q_cellnw(:,:,:)
53  REAL*8, ALLOCATABLE :: q_cellne(:,:,:)
55  REAL*8, ALLOCATABLE :: q_cellsw(:,:,:)
57  REAL*8, ALLOCATABLE :: q_cellse(:,:,:)
58 
60  REAL*8, ALLOCATABLE :: q1max(:,:)
61 
63  REAL*8, ALLOCATABLE :: a_interface_xneg(:,:,:)
65  REAL*8, ALLOCATABLE :: a_interface_xpos(:,:,:)
67  REAL*8, ALLOCATABLE :: a_interface_yneg(:,:,:)
69  REAL*8, ALLOCATABLE :: a_interface_ypos(:,:,:)
71  REAL*8, ALLOCATABLE :: h_interface_x(:,:,:)
73  REAL*8, ALLOCATABLE :: h_interface_y(:,:,:)
75  REAL*8, ALLOCATABLE :: qp(:,:,:)
76 
78  REAL*8, ALLOCATABLE :: source_xy(:,:)
79 
80  LOGICAL, ALLOCATABLE :: solve_mask(:,:) , solve_mask0(:,:)
81 
83  REAL*8 :: dt
84 
85  LOGICAL, ALLOCATABLE :: mask22(:,:) , mask21(:,:) , mask11(:,:) , mask12(:,:)
86 
87  INTEGER :: i_rk
88 
90  REAL*8, ALLOCATABLE :: a_tilde_ij(:,:)
92  REAL*8, ALLOCATABLE :: a_dirk_ij(:,:)
93 
95  REAL*8, ALLOCATABLE :: omega_tilde(:)
96 
98  REAL*8, ALLOCATABLE :: omega(:)
99 
101  REAL*8, ALLOCATABLE :: a_tilde(:)
102 
104  REAL*8, ALLOCATABLE :: a_dirk(:)
105 
107  REAL*8 :: a_diag
108 
110  REAL*8, ALLOCATABLE :: q_rk(:,:,:,:)
111 
113  REAL*8, ALLOCATABLE :: divflux(:,:,:,:)
114 
116  REAL*8, ALLOCATABLE :: nh(:,:,:,:)
117 
119  REAL*8, ALLOCATABLE :: si_nh(:,:,:,:)
120 
122  REAL*8, ALLOCATABLE :: expl_terms(:,:,:,:)
123 
125  REAL*8, ALLOCATABLE :: divfluxj(:,:)
126 
128  REAL*8, ALLOCATABLE :: nhj(:,:)
129 
131  REAL*8, ALLOCATABLE :: expl_terms_j(:,:)
132 
134  REAL*8, ALLOCATABLE :: si_nhj(:,:)
135 
137  LOGICAL :: normalize_q
138 
140  LOGICAL :: normalize_f
141 
143  LOGICAL :: opt_search_nl
144 
146  REAL*8, ALLOCATABLE :: residual_term(:,:,:)
147 
148 
149  REAL*8 :: t_imex1,t_imex2
150 
151 CONTAINS
152 
153  !******************************************************************************
155  !
158  !
162  !
163  !******************************************************************************
164 
165  SUBROUTINE allocate_solver_variables
166 
167  IMPLICIT NONE
168 
169  REAL*8 :: gamma , delta
170 
171  INTEGER :: i,j
172 
173  ALLOCATE( q( n_vars , comp_cells_x , comp_cells_y ) , q0( n_vars , &
174  comp_cells_x , comp_cells_y ) )
175 
176  ALLOCATE( q1max( comp_cells_x , comp_cells_y ) )
177 
178  ALLOCATE( q_fv( n_vars , comp_cells_x , comp_cells_y ) )
179 
180  ALLOCATE( q_interfacel( n_vars , comp_interfaces_x, comp_cells_y ) )
181  ALLOCATE( q_interfacer( n_vars , comp_interfaces_x, comp_cells_y ) )
182  ALLOCATE( a_interface_xneg( n_eqns , comp_interfaces_x, comp_cells_y ) )
183  ALLOCATE( a_interface_xpos( n_eqns , comp_interfaces_x, comp_cells_y ) )
184  ALLOCATE( h_interface_x( n_eqns , comp_interfaces_x, comp_cells_y ) )
185 
186 
187  ALLOCATE( q_interfaceb( n_vars , comp_cells_x, comp_interfaces_y ) )
188  ALLOCATE( q_interfacet( n_vars , comp_cells_x, comp_interfaces_y ) )
189  ALLOCATE( a_interface_yneg( n_eqns , comp_cells_x, comp_interfaces_y ) )
190  ALLOCATE( a_interface_ypos( n_eqns , comp_cells_x, comp_interfaces_y ) )
191 
192  ALLOCATE( q_cellnw( n_vars , comp_cells_x , comp_cells_y ) )
193  ALLOCATE( q_cellne( n_vars , comp_cells_x , comp_cells_y ) )
194  ALLOCATE( q_cellsw( n_vars , comp_cells_x , comp_cells_y ) )
195  ALLOCATE( q_cellse( n_vars , comp_cells_x , comp_cells_y ) )
196 
197  ALLOCATE( h_interface_y( n_eqns , comp_cells_x, comp_interfaces_y ) )
198 
199  ALLOCATE( solve_mask( comp_cells_x , comp_cells_y ) )
200  ALLOCATE( solve_mask0( comp_cells_x , comp_cells_y ) )
201 
202  ALLOCATE( qp( n_vars , comp_cells_x , comp_cells_y ) )
203 
204  ALLOCATE( source_xy( comp_cells_x , comp_cells_y ) )
205 
206 
207  ALLOCATE( a_tilde_ij(n_rk,n_rk) )
208  ALLOCATE( a_dirk_ij(n_rk,n_rk) )
209  ALLOCATE( omega_tilde(n_rk) )
210  ALLOCATE( omega(n_rk) )
211 
212 
213  ! Allocate the logical arrays defining the implicit parts of the system
214  ALLOCATE( mask22(n_eqns,n_eqns) )
215  ALLOCATE( mask21(n_eqns,n_eqns) )
216  ALLOCATE( mask11(n_eqns,n_eqns) )
217  ALLOCATE( mask12(n_eqns,n_eqns) )
218 
219  ! Initialize the logical arrays with all false (everything is implicit)
220  mask11(1:n_eqns,1:n_eqns) = .false.
221  mask12(1:n_eqns,1:n_eqns) = .false.
222  mask22(1:n_eqns,1:n_eqns) = .false.
223  mask21(1:n_eqns,1:n_eqns) = .false.
224 
225  ! Set to .TRUE. the elements not corresponding to equations and variables to
226  ! be solved implicitly
227  DO i = 1,n_eqns
228 
229  DO j = 1,n_eqns
230 
231  IF ( .NOT.implicit_flag(i) .AND. .NOT.implicit_flag(j) ) &
232  mask11(j,i) = .true.
233  IF ( implicit_flag(i) .AND. .NOT.implicit_flag(j) ) &
234  mask12(j,i) = .true.
235  IF ( implicit_flag(i) .AND. implicit_flag(j) ) &
236  mask22(j,i) = .true.
237  IF ( .NOT.implicit_flag(i) .AND. implicit_flag(j) ) &
238  mask21(j,i) = .true.
239 
240  END DO
241 
242  END DO
243 
244  ! Initialize the coefficients for the IMEX Runge-Kutta scheme
245  ! Please note that with respect to the schemes described in Pareschi & Russo
246  ! (2000) we do not have the coefficient vectors c_tilde and c, because the
247  ! explicit and implicit terms do not depend explicitly on time.
248 
249  ! Explicit part coefficients (a_tilde_ij=0 for j>=i)
250  a_tilde_ij = 0.d0
251 
252  ! Weight coefficients of the explicit part in the final assemblage
253  omega_tilde = 0.d0
254 
255  ! Implicit part coefficients (a_dirk_ij=0 for j>i)
256  a_dirk_ij = 0.d0
257 
258  ! Weight coefficients of the explicit part in the final assemblage
259  omega = 0.d0
260 
261  gamma = 1.d0 - 1.d0 / sqrt(2.d0)
262  delta = 1.d0 - 1.d0 / ( 2.d0 * gamma )
263 
264  IF ( n_rk .EQ. 1 ) THEN
265 
266  a_tilde_ij(1,1) = 1.d0
267 
268  omega_tilde(1) = 1.d0
269 
270  a_dirk_ij(1,1) = 0.d0
271 
272  omega(1) = 0.d0
273 
274  ELSEIF ( n_rk .EQ. 2 ) THEN
275 
276  a_tilde_ij(2,1) = 1.0d0
277 
278  omega_tilde(1) = 1.0d0
279  omega_tilde(2) = 0.0d0
280 
281  a_dirk_ij(2,2) = 1.0d0
282 
283  omega(1) = 0.d0
284  omega(2) = 1.d0
285 
286 
287  ELSEIF ( n_rk .EQ. 3 ) THEN
288 
289  ! Tableau for the IMEX-SSP(3,3,2) Stiffly Accurate Scheme
290  ! from Pareschi & Russo (2005), Table IV
291 
292  a_tilde_ij(2,1) = 0.5d0
293  a_tilde_ij(3,1) = 0.5d0
294  a_tilde_ij(3,2) = 0.5d0
295 
296  omega_tilde(1) = 1.0d0 / 3.0d0
297  omega_tilde(2) = 1.0d0 / 3.0d0
298  omega_tilde(3) = 1.0d0 / 3.0d0
299 
300  a_dirk_ij(1,1) = 0.25d0
301  a_dirk_ij(2,2) = 0.25d0
302  a_dirk_ij(3,1) = 1.0d0 / 3.0d0
303  a_dirk_ij(3,2) = 1.0d0 / 3.0d0
304  a_dirk_ij(3,3) = 1.0d0 / 3.0d0
305 
306  omega(1) = 1.0d0 / 3.0d0
307  omega(2) = 1.0d0 / 3.0d0
308  omega(3) = 1.0d0 / 3.0d0
309 
310 
311 
312  ELSEIF ( n_rk .EQ. 4 ) THEN
313 
314  ! LRR(3,2,2) from Table 3 in Pareschi & Russo (2000)
315 
316  a_tilde_ij(2,1) = 0.5d0
317  a_tilde_ij(3,1) = 1.d0 / 3.d0
318  a_tilde_ij(4,2) = 1.0d0
319 
320  omega_tilde(1) = 0.d0
321  omega_tilde(2) = 1.0d0
322  omega_tilde(3) = 0.0d0
323  omega_tilde(4) = 0.d0
324 
325  a_dirk_ij(2,2) = 0.5d0
326  a_dirk_ij(3,3) = 1.0d0 / 3.0d0
327  a_dirk_ij(4,3) = 0.75d0
328  a_dirk_ij(4,4) = 0.25d0
329 
330  omega(1) = 0.d0
331  omega(2) = 0.d0
332  omega(3) = 0.75d0
333  omega(4) = 0.25d0
334 
335  END IF
336 
337  ALLOCATE( a_tilde(n_rk) )
338  ALLOCATE( a_dirk(n_rk) )
339 
340  ALLOCATE( q_rk( n_vars , comp_cells_x , comp_cells_y , n_rk ) )
341  ALLOCATE( divflux( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )
342  ALLOCATE( nh( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )
343  ALLOCATE( si_nh( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )
344 
345  ALLOCATE( expl_terms( n_eqns , comp_cells_x , comp_cells_y , n_rk ) )
346 
347  ALLOCATE( divfluxj(n_eqns,n_rk) )
348  ALLOCATE( nhj(n_eqns,n_rk) )
349  ALLOCATE( si_nhj(n_eqns,n_rk) )
350 
351  ALLOCATE( expl_terms_j(n_eqns,n_rk) )
352 
353  ALLOCATE( residual_term( n_vars , comp_cells_x , comp_cells_y ) )
354 
355  END SUBROUTINE allocate_solver_variables
356 
357  !******************************************************************************
359  !
362  !
366  !
367  !******************************************************************************
368 
369  SUBROUTINE deallocate_solver_variables
370 
371  DEALLOCATE( q , q0 )
372 
373  DEALLOCATE( q1max )
374 
375  DEALLOCATE( q_fv )
376 
377  DEALLOCATE( q_interfacel )
378  DEALLOCATE( q_interfacer )
379  DEALLOCATE( q_interfaceb )
380  DEALLOCATE( q_interfacet )
381 
382  DEALLOCATE( q_cellnw )
383  DEALLOCATE( q_cellne )
384  DEALLOCATE( q_cellsw )
385  DEALLOCATE( q_cellse )
386 
387  DEALLOCATE( a_interface_xneg )
388  DEALLOCATE( a_interface_xpos )
389  DEALLOCATE( a_interface_yneg )
390  DEALLOCATE( a_interface_ypos )
391 
392  DEALLOCATE( h_interface_x )
393  DEALLOCATE( h_interface_y )
394 
395  DEALLOCATE( solve_mask , solve_mask0 )
396 
397  Deallocate( qp )
398 
399  DEALLOCATE( source_xy )
400 
401  DEALLOCATE( a_tilde_ij )
402  DEALLOCATE( a_dirk_ij )
403  DEALLOCATE( omega_tilde )
404  DEALLOCATE( omega )
405 
406  DEALLOCATE( implicit_flag )
407 
408  DEALLOCATE( a_tilde )
409  DEALLOCATE( a_dirk )
410 
411  DEALLOCATE( q_rk )
412  DEALLOCATE( divflux )
413  DEALLOCATE( nh )
414  DEALLOCATE( si_nh )
415  DEALLOCATE( expl_terms )
416 
417  DEALLOCATE( divfluxj )
418  DEALLOCATE( nhj )
419  DEALLOCATE( si_nhj )
420  DEALLOCATE( expl_terms_j )
421 
422  DEALLOCATE( mask22 , mask21 , mask11 , mask12 )
423 
424  DEALLOCATE( residual_term )
425 
426  END SUBROUTINE deallocate_solver_variables
427 
428 
429  !******************************************************************************
431  !
435  !
439  !
440  !******************************************************************************
441 
442  SUBROUTINE check_solve
443 
444  IMPLICIT NONE
445 
446  INTEGER :: i
447 
448  solve_mask0(1:comp_cells_x,1:comp_cells_y) = .false.
449 
450  WHERE ( q(1,:,:) .GT. 0.d0 ) solve_mask0 = .true.
451 
452  solve_mask = solve_mask0
453 
454  DO i = 1,n_rk
455 
456  solve_mask(1+i:comp_cells_x,:) = solve_mask(1+i:comp_cells_x,:) .OR. &
457  solve_mask(1:comp_cells_x-i,:)
458 
459  solve_mask(1:comp_cells_x-i,:) = solve_mask(1:comp_cells_x-i,:) .OR. &
460  solve_mask(1+i:comp_cells_x,:)
461 
462  solve_mask(:,1+i:comp_cells_y) = solve_mask(:,1+i:comp_cells_y) .OR. &
463  solve_mask(:,1:comp_cells_y-i)
464 
465  solve_mask(:,1:comp_cells_y-i) = solve_mask(:,1:comp_cells_y-i) .OR. &
466  solve_mask(:,1+i:comp_cells_y)
467 
468  END DO
469 
470 
471  END SUBROUTINE check_solve
472 
473  !******************************************************************************
475  !
479  !
483  !
484  !******************************************************************************
485 
486  SUBROUTINE timestep
487 
488  ! External variables
489  USE geometry_2d, ONLY : dx,dy
490  USE parameters_2d, ONLY : max_dt , cfl
491 
492  ! External procedures
493  USE constitutive_2d, ONLY : eval_local_speeds_x, eval_local_speeds_y
494 
495  IMPLICIT none
496 
497  REAL*8 :: vel_max(n_vars)
498  REAL*8 :: vel_min(n_vars)
499  REAL*8 :: vel_j
500  REAL*8 :: dt_cfl
501  REAL*8 :: qj(n_vars)
502 
503  INTEGER :: j,k
504 
505  dt = max_dt
506 
507  IF ( cfl .NE. -1.d0 ) THEN
508 
509  DO j = 1,comp_cells_x
510 
511  DO k = 1,comp_cells_y
512 
513  qj = q( 1:n_vars , j , k )
514 
515  IF ( comp_cells_x .GT. 1 ) THEN
516 
517  ! x direction
518  CALL eval_local_speeds_x( qj , vel_min , vel_max )
519 
520  vel_j = max( maxval(abs(vel_min)) , maxval(abs(vel_max)) )
521 
522  dt_cfl = cfl * dx / vel_j
523 
524  dt = min( dt , dt_cfl )
525 
526  END IF
527 
528  IF ( comp_cells_y .GT. 1 ) THEN
529 
530  ! y direction
531  CALL eval_local_speeds_y( qj , vel_min , vel_max )
532 
533  vel_j = max( maxval(abs(vel_min)) , maxval(abs(vel_max)) )
534 
535  dt_cfl = cfl * dy / vel_j
536 
537  dt = min( dt , dt_cfl )
538 
539  END IF
540 
541  ENDDO
542 
543  END DO
544 
545  END IF
546 
547  END SUBROUTINE timestep
548 
549 
550  !*****************************************************************************
552  !
556  !
560  !
561  !*****************************************************************************
562 
563  SUBROUTINE timestep2
564 
565  ! External variables
566  USE geometry_2d, ONLY : dx,dy
567  USE parameters_2d, ONLY : max_dt , cfl
568 
569  IMPLICIT none
570 
571  REAL*8 :: dt_cfl
572 
573  REAL*8 :: a_interface_x_max(n_eqns,comp_interfaces_x,comp_cells_y)
574  REAL*8 :: a_interface_y_max(n_eqns,comp_cells_x,comp_interfaces_y)
575  REAL*8 :: dt_interface_x, dt_interface_y
576 
577  INTEGER :: i,j,k
578 
579  dt = max_dt
580 
581  IF ( cfl .NE. -1.d0 ) THEN
582 
583  CALL reconstruction
584 
585  CALL eval_speeds
586 
587  DO i=1,n_vars
588 
589  a_interface_x_max(i,:,:) = &
590  max( a_interface_xpos(i,:,:) , -a_interface_xneg(i,:,:) )
591 
592  a_interface_y_max(i,:,:) = &
593  max( a_interface_ypos(i,:,:) , -a_interface_yneg(i,:,:) )
594 
595  END DO
596 
597  DO j = 1,comp_cells_x
598 
599  DO k = 1,comp_cells_y
600 
601  dt_interface_x = cfl * dx / max( maxval(a_interface_x_max(:,j,k)) &
602  , maxval(a_interface_x_max(:,j+1,k)) )
603 
604  dt_interface_y = cfl * dy / max( maxval(a_interface_y_max(:,j,k)) &
605  , maxval(a_interface_y_max(:,j,k+1)) )
606 
607  dt_cfl = min( dt_interface_x , dt_interface_y )
608 
609  dt = min(dt,dt_cfl)
610 
611  ENDDO
612 
613  END DO
614 
615  END IF
616 
617  END SUBROUTINE timestep2
618 
619  !******************************************************************************
621  !
626  !
630  !
631  !******************************************************************************
632 
633  SUBROUTINE imex_rk_solver
634 
635  USE constitutive_2d, ONLY : eval_nonhyperbolic_terms
636 
637  USE constitutive_2d, ONLY : eval_nh_semi_impl_terms
638 
639  USE constitutive_2d, ONLY : qc_to_qp
640 
641  IMPLICIT NONE
642 
643  REAL*8 :: q_si(n_vars)
644  REAL*8 :: q_guess(n_vars)
645  INTEGER :: j,k
646 
647  REAL*8 :: h_new
648 
649  ! Initialization of the solution guess
650  q0( 1:n_vars , 1:comp_cells_x , 1:comp_cells_y ) = &
651  q( 1:n_vars , 1:comp_cells_x , 1:comp_cells_y )
652 
653  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solver, imex_RK_solver: beginning'
654 
655  ! Initialization of the variables for the Runge-Kutta scheme
656  q_rk(1:n_vars,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0
657 
658  divflux(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0
659 
660  nh(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0
661 
662  si_nh(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0
663 
664  expl_terms(1:n_eqns,1:comp_cells_x,1:comp_cells_y,1:n_rk) = 0.d0
665 
666  runge_kutta:DO i_rk = 1,n_rk
667 
668  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solver, imex_RK_solver: i_RK',i_rk
669 
670  ! define the explicits coefficients for the i-th step of the Runge-Kutta
671  a_tilde = 0.d0
672  a_dirk = 0.d0
673 
674  ! in the first step of the RK scheme all the coefficients remain to 0
675  a_tilde(1:i_rk-1) = a_tilde_ij(i_rk,1:i_rk-1)
676  a_dirk(1:i_rk-1) = a_dirk_ij(i_rk,1:i_rk-1)
677 
678  ! define the implicit coefficient for the i-th step of the Runge-Kutta
679  a_diag = a_dirk_ij(i_rk,i_rk)
680 
681  CALL cpu_time(t_imex1)
682 
683  loop_over_ycells:DO k = 1,comp_cells_y
684 
685  loop_over_xcells:DO j = 1,comp_cells_x
686 
687  IF ( verbose_level .GE. 2 ) THEN
688 
689  WRITE(*,*) 'solver, imex_RK_solver: j',j,k
690 
691  END IF
692 
693  ! initialize the RK step
694  IF ( i_rk .EQ. 1 ) THEN
695 
696  ! solution from the previous time step
697  q_guess(1:n_vars) = q0( 1:n_vars , j , k)
698 
699  ELSE
700 
701  ! solution from the previous RK step
702  !q_guess(1:n_vars) = q_rk( 1:n_vars , j , k , MAX(1,i_RK-1) )
703 
704  END IF
705 
706  ! RK explicit hyperbolic terms for the volume (i,j)
707  divfluxj(1:n_eqns,1:n_rk) = divflux( 1:n_eqns , j , k , 1:n_rk )
708 
709  ! RK implicit terms for the volume (i,j)
710  nhj(1:n_eqns,1:n_rk) = nh( 1:n_eqns , j , k , 1:n_rk )
711 
712  ! RK semi-implicit terms for the volume (i,j)
713  si_nhj(1:n_eqns,1:n_rk) = si_nh( 1:n_eqns , j , k , 1:n_rk )
714 
715  ! RK additional explicit terms for the volume (i,j)
716  expl_terms_j(1:n_eqns,1:n_rk) = expl_terms( 1:n_eqns,j,k,1:n_rk )
717 
718  ! New solution at the i_RK step without the implicit and
719  ! semi-implicit term
720  q_fv( 1:n_vars , j , k ) = q0( 1:n_vars , j , k ) &
721  - dt * (matmul( divfluxj(1:n_eqns,1:i_rk) &
722  + expl_terms_j(1:n_eqns,1:i_rk) , a_tilde(1:i_rk) ) &
723  - matmul( nhj(1:n_eqns,1:i_rk) + si_nhj(1:n_eqns,1:i_rk) , &
724  a_dirk(1:i_rk) ) )
725 
726  IF ( verbose_level .GE. 2 ) THEN
727 
728  WRITE(*,*) 'q_guess',q_guess
729  CALL qc_to_qp( q_guess , b_cent(j,k) , qp(1:n_vars,j,k) )
730  WRITE(*,*) 'q_guess: qp',qp(1:n_vars,j,k)
731 
732  END IF
733 
734  adiag_pos:IF ( a_diag .NE. 0.d0 ) THEN
735 
736  pos_thick:IF ( q_fv(1,j,k) .GT. 0.d0 ) THEN
737 
738  ! WRITE(*,*) 'SOLVER_2D, imex_RK_solver: q_fv', q_fv(1:n_vars,j,k )
739 
740  ! Eval the semi-implicit discontinuous terms
741  CALL eval_nh_semi_impl_terms( grav_surf(j,k) , &
742  r_qj = q_fv( 1:n_vars , j , k ) , &
743  r_nh_semi_impl_term = si_nh(1:n_eqns,j,k,i_rk) )
744 
745  ! WRITE(*,*) 'SOLVER_2D, imex_RK_solver: SI_NH', SI_NH(1:n_eqns,j,k,i_RK)
746 
747  si_nhj(1:n_eqns,i_rk) = si_nh( 1:n_eqns,j,k,i_rk )
748 
749  ! Assemble the initial guess for the implicit solver
750  q_si(1:n_vars) = q_fv(1:n_vars,j,k ) + dt * a_diag * &
751  si_nh(1:n_eqns,j,k,i_rk)
752 
753  ! WRITE(*,*) 'SOLVER_2D, imex_RK_solver: q_si', q_si(1:n_vars )
754 
755  IF ( q_fv(2,j,k)**2 + q_fv(3,j,k)**2 .EQ. 0.d0 ) THEN
756 
757  !Case 1: if the velocity was null, then it must stay null
758  q_si(2:3) = 0.d0
759 
760  ELSEIF ( ( q_si(2)*q_fv(2,j,k) .LT. 0.d0 ) .OR. &
761  ( q_si(3)*q_fv(3,j,k) .LT. 0.d0 ) ) THEN
762 
763  ! If the semi-impl. friction term changed the sign of the
764  ! velocity then set it to zero
765  q_si(2:3) = 0.d0
766 
767  ELSE
768 
769  ! Align the velocity vector with previous one
770  q_si(2:3) = dsqrt( q_si(2)**2 + q_si(3)**2 ) * &
771  q_fv(2:3,j,k) / dsqrt( q_fv(2,j,k)**2 &
772  + q_fv(3,j,k)**2 )
773 
774  END IF
775 
776  ! Update the semi-implicit term accordingly with the
777  ! corrections above
778  si_nh(1:n_eqns,j,k,i_rk) = ( q_si(1:n_vars) - &
779  q_fv(1:n_vars,j,k ) ) / ( dt*a_diag )
780 
781  si_nhj(1:n_eqns,i_rk) = si_nh( 1:n_eqns,j,k,i_rk )
782 
783  ! Initialize the guess for the NR solver
784  q_guess(1:n_vars) = q_si(1:n_vars)
785 
786  ! WRITE(*,*) 'SOLVER_2D, imex_RK_solver: q_guess', q_guess(1:n_vars )
787 
788  ! WRITE(*,*) 'SOLVER_2D, imex_RK_solver: j,k',j,k
789 
790  ! Solve the implicit system to find the solution at the
791  ! i_RK step of the IMEX RK procedure
792  CALL solve_rk_step( b_cent(j,k) , q_guess(1:n_vars) , &
793  q0(1:n_vars,j,k ) , a_tilde , a_dirk , a_diag )
794 
795  IF ( comp_cells_y .EQ. 1 ) THEN
796 
797  q_guess(3) = 0.d0
798 
799  END IF
800 
801  IF ( comp_cells_x .EQ. 1 ) THEN
802 
803  q_guess(2) = 0.d0
804 
805  END IF
806 
807  ! Eval and store the implicit term at the i_RK step
808  CALL eval_nonhyperbolic_terms( r_qj =q_guess , &
809  r_nh_term_impl = nh(1:n_eqns,j,k,i_rk) )
810 
811  IF ( q_si(2)**2 + q_si(3)**2 .EQ. 0.d0 ) THEN
812 
813  q_guess(2:3) = 0.d0
814 
815  ELSEIF ( ( q_guess(2)*q_si(2) .LE. 0.d0 ) .AND. &
816  ( q_guess(3)*q_si(3) .LE. 0.d0 ) ) THEN
817 
818  ! If the impl. friction term changed the sign of the
819  ! velocity then set it to zero
820  q_guess(2:3) = 0.d0
821 
822  ELSE
823 
824  ! Align the velocity vector with previous one
825  q_guess(2:3) = dsqrt( q_guess(2)**2 + q_guess(3)**2 ) * &
826  q_si(2:3) / dsqrt( q_si(2)**2 + q_si(3)**2 )
827 
828  END IF
829 
830  ELSE
831 
832  ! If h=0 nothing has to be changed
833  q_guess(1:n_vars) = q_fv( 1:n_vars , j , k )
834  q_si(1:n_vars) = q_fv( 1:n_vars , j , k )
835  si_nh(1:n_eqns,j,k,i_rk) = 0.d0
836  nh(1:n_eqns,j,k,i_rk) = 0.d0
837 
838  END IF pos_thick
839 
840  END IF adiag_pos
841 
842  ! Check the sign of the flow thickness
843  h_new = q_guess(1)
844 
845  IF ( j .EQ. 0 ) THEN
846 
847  WRITE(*,*) 'h',q_guess(1)
848  READ(*,*)
849 
850  END IF
851 
852  ! IF ( h_new .LT. 1.D-10 ) THEN
853  !
854  ! q_guess(1) = B_cent(j,k)
855  ! q_guess(2:n_vars) = 0.D0
856  !
857  ! END IF
858 
859  IF ( a_diag .NE. 0.d0 ) THEN
860 
861  ! Update the viscous term with the correction on the new velocity
862  nh(1:n_vars,j,k,i_rk) = ( q_guess(1:n_vars) - q_si(1:n_vars)) &
863  / ( dt*a_diag )
864 
865  END IF
866 
867  ! Store the solution at the end of the i_RK step
868  q_rk( 1:n_vars , j , k , i_rk ) = q_guess
869 
870  IF ( verbose_level .GE. 2 ) THEN
871 
872  WRITE(*,*) 'imex_RK_solver: qc',q_guess
873  CALL qc_to_qp( q_guess, b_cent(j,k) , qp(1:n_vars,j,k) )
874  WRITE(*,*) 'imex_RK_solver: qp',qp(1:n_vars,j,k)
875  READ(*,*)
876 
877  END IF
878 
879  END DO loop_over_xcells
880 
881  ENDDO loop_over_ycells
882 
883  CALL cpu_time(t_imex2)
884  ! WRITE(*,*) 'Time taken by implicit',t_imex2-t_imex1,'seconds'
885 
886  IF ( omega_tilde(i_rk) .GT. 0.d0 ) THEN
887 
888  ! Eval and store the explicit hyperbolic (fluxes) terms
889  CALL eval_hyperbolic_terms( q_rk(1:n_vars,1:comp_cells_x,1:comp_cells_y, &
890  i_rk) , divflux(1:n_eqns,1:comp_cells_x,1:comp_cells_y,i_rk) )
891 
892  CALL cpu_time(t_imex1)
893  ! WRITE(*,*) 'Time taken by explicit',t_imex1-t_imex2,'seconds'
894 
895  ! Eval and store the other explicit terms (e.g. gravity or viscous forces)
896  CALL eval_explicit_terms( q_rk(1:n_vars,1:comp_cells_x,1:comp_cells_y, &
897  i_rk) , expl_terms(1:n_eqns,1:comp_cells_x,1:comp_cells_y,i_rk) )
898 
899  CALL cpu_time(t_imex1)
900  ! WRITE(*,*) 'Time taken by explicit',t_imex1-t_imex2,'seconds'
901 
902  END IF
903 
904  IF ( verbose_level .GE. 1 ) THEN
905 
906  WRITE(*,*) 'div_flux(2),div_flux(3),expl_terms(2),expl_terms(3)'
907 
908  DO k = 1,comp_cells_y
909 
910  DO j = 1,comp_cells_x
911 
912  WRITE(*,*) divflux(2,j,k,i_rk) , divflux(3,j,k,i_rk) , &
913  expl_terms(2,j,k,i_rk) , expl_terms(3,j,k,i_rk)
914 
915  ENDDO
916 
917  END DO
918 
919  READ(*,*)
920 
921  END IF
922 
923  END DO runge_kutta
924 
925  DO k = 1,comp_cells_y
926 
927  DO j = 1,comp_cells_x
928 
929  residual_term(1:n_vars,j,k) = matmul( divflux(1:n_eqns,j,k,1:n_rk) &
930  + expl_terms(1:n_eqns,j,k,1:n_rk) , omega_tilde ) - &
931  matmul( nh(1:n_eqns,j,k,1:n_rk) + si_nh(1:n_eqns,j,k,1:n_rk) ,&
932  omega )
933 
934  ENDDO
935 
936  END DO
937 
938  assemble_sol_loop_x:DO k = 1,comp_cells_y
939 
940  assemble_sol_loop_y:DO j = 1,comp_cells_x
941 
942  IF ( verbose_level .GE. 1 ) THEN
943 
944  WRITE(*,*) 'cell jk =',j,k
945  WRITE(*,*) 'before imex_RK_solver: qc',q0(1:n_vars,j,k)
946  CALL qc_to_qp(q0(1:n_vars,j,k) , b_cent(j,k) , qp(1:n_vars,j,k))
947  WRITE(*,*) 'before imex_RK_solver: qp',qp(1:n_vars,j,k)
948 
949  END IF
950 
951  IF ( ( sum(abs( omega_tilde(:)-a_tilde_ij(n_rk,:))) .EQ. 0.d0 ) &
952  .AND. ( sum(abs(omega(:)-a_dirk_ij(n_rk,:))) .EQ. 0.d0 ) ) THEN
953 
954  ! The assembling coeffs are equal to the last step of the RK scheme
955  q(1:n_vars,j,k) = q_rk(1:n_vars,j,k,n_rk)
956 
957  ELSE
958 
959  ! The assembling coeffs are different
960  q(1:n_vars,j,k) = q0(1:n_vars,j,k) - dt*residual_term(1:n_vars,j,k)
961 
962  END IF
963 
964  IF ( ( j .EQ. -1 ) .AND. ( k .EQ. 36 ) ) THEN
965 
966  WRITE(*,*) 'after assemble new solution'
967  WRITE(*,*) 'j,k',j,k
968  WRITE(*,*) 'q_old(1,j,k),q_new(1,j,k)',q0(1,j,k), q(1,j,k)
969  WRITE(*,*) 'q_old(4,j,k),q_new(4,j,k)',q0(4,j,k), q(4,j,k)
970  WRITE(*,*) 'divFlux(1,j,k,1:n_RK)',divflux(1,j,k,1:n_rk)
971  WRITE(*,*) 'dt',dt
972 
973  END IF
974 
975  negative_thickness_check:IF ( q(1,j,k) .LT. 0.d0 ) THEN
976 
977  IF ( q(1,j,k) .GT. -1.d-7 ) THEN
978 
979  q(1,j,k) = 0.d0
980  q(2:n_vars,j,k) = 0.d0
981 
982  ELSE
983 
984  WRITE(*,*) 'j,k,n_RK',j,k,n_rk
985  WRITE(*,*) 'dt',dt
986 
987  WRITE(*,*) 'before imex_RK_solver: qc',q0(1:n_vars,j,k)
988  CALL qc_to_qp(q0(1:n_vars,j,k) , b_cent(j,k) , qp(1:n_vars,j,k))
989  WRITE(*,*) 'before imex_RK_solver: qp',qp(1:n_vars,j,k)
990  WRITE(*,*) 'h old',q0(1,j,k)
991  WRITE(*,*) 'h new',q(1,j,k)
992  WRITE(*,*) 'B_cent(j,k)',b_cent(j,k)
993  WRITE(*,*) 'B_stag_x(j:j+1,k)',b_stag_x(j:j+1,k)
994  WRITE(*,*) 'B_stag_y(j,k:k+1)',b_stag_y(j,k:k+1)
995 
996  WRITE(*,*) 'hS',q_interfacet(1,j,k)
997  WRITE(*,*) 'hE',q_interfacer(1,j,k)
998 
999  READ(*,*)
1000 
1001  END IF
1002 
1003  END IF negative_thickness_check
1004 
1005  IF ( solid_transport_flag ) THEN
1006 
1007  IF ( q(4,j,k) .GT. q(1,j,k) ) THEN
1008 
1009  IF ( q(4,j,k)-q(1,j,k) .LT. 1.d-10 ) THEN
1010 
1011  q(4,j,k) = q(1,j,k)
1012 
1013  ELSE
1014 
1015  WRITE(*,*) 'j,k,n_RK',j,k,n_rk
1016  WRITE(*,*) 'dt',dt
1017  WRITE(*,*) ' B_cent(j,k)', b_cent(j,k)
1018 
1019  WRITE(*,*) 'before imex_RK_solver: qc',q0(1:n_vars,j,k)
1020  CALL qc_to_qp(q0(1:n_vars,j,k) , b_cent(j,k) , qp(1:n_vars,j,k))
1021  WRITE(*,*) 'before imex_RK_solver: qp',qp(1:n_vars,j,k)
1022 
1023  CALL qc_to_qp(q(1:n_vars,j,k) , b_cent(j,k) , qp(1:n_vars,j,k))
1024 
1025  WRITE(*,*) 'after imex_RK_solver: qc',q(1:n_vars,j,k)
1026  WRITE(*,*) 'after imex_RK_solver: qp',qp(1:n_vars,j,k)
1027 
1028 
1029  WRITE(*,*) 'h old',q0(1,j,k)
1030  WRITE(*,*) 'h new',q(1,j,k)
1031  WRITE(*,*) 'alphas old', q0(4,j,k) / q0(1,j,k)
1032  WRITE(*,*) 'alphas new', q(4,j,k) / q(1,j,k)
1033  READ(*,*)
1034 
1035  END IF
1036 
1037  END IF
1038 
1039  IF ( verbose_level .GE. 1 ) THEN
1040 
1041  WRITE(*,*) 'h new',q(1,j,k)
1042  READ(*,*)
1043 
1044  END IF
1045 
1046  END IF
1047 
1048  ENDDO assemble_sol_loop_y
1049 
1050  END DO assemble_sol_loop_x
1051 
1052  RETURN
1053 
1054  END SUBROUTINE imex_rk_solver
1055 
1056  !******************************************************************************
1058  !
1065  !
1073  !
1077  !
1078  !******************************************************************************
1079 
1080  SUBROUTINE solve_rk_step( Bj, qj, qj_old, a_tilde , a_dirk , a_diag )
1081 
1082  USE parameters_2d, ONLY : max_nl_iter , tol_rel , tol_abs
1083 
1084  USE constitutive_2d, ONLY : qc_to_qp
1085 
1086  IMPLICIT NONE
1087 
1088  REAL*8, INTENT(IN) :: Bj
1089  REAL*8, INTENT(INOUT) :: qj(n_vars)
1090  REAL*8, INTENT(IN) :: qj_old(n_vars)
1091  REAL*8, INTENT(IN) :: a_tilde(n_rk)
1092  REAL*8, INTENT(IN) :: a_dirk(n_rk)
1093  REAL*8, INTENT(IN) :: a_diag
1094 
1095  REAL*8 :: qj_init(n_vars)
1096 
1097  REAL*8 :: qj_org(n_vars) , qj_rel(n_vars)
1098 
1099  REAL*8 :: left_matrix(n_eqns,n_vars)
1100  REAL*8 :: right_term(n_eqns)
1101 
1102  REAL*8 :: scal_f
1103 
1104  REAL*8 :: coeff_f(n_eqns)
1105 
1106  REAL*8 :: qj_rel_NR_old(n_vars)
1107  REAL*8 :: scal_f_old
1108  REAL*8 :: desc_dir(n_vars)
1109  REAL*8 :: grad_f(n_vars)
1110  REAL*8 :: mod_desc_dir
1111 
1112  INTEGER :: pivot(n_vars)
1113 
1114  REAL*8 :: left_matrix_small22(n_nh,n_nh)
1115  REAL*8 :: left_matrix_small21(n_eqns-n_nh,n_nh)
1116  REAL*8 :: left_matrix_small11(n_eqns-n_nh,n_vars-n_nh)
1117  REAL*8 :: left_matrix_small12(n_nh,n_vars-n_nh)
1118 
1119  REAL*8 :: desc_dir_small2(n_nh)
1120  INTEGER :: pivot_small2(n_nh)
1121 
1122  REAL*8 :: desc_dir_small1(n_vars-n_nh)
1123 
1124  INTEGER :: ok
1125 
1126  INTEGER :: i
1127  INTEGER :: nl_iter
1128 
1129  REAL*8, PARAMETER :: STPMX=100.d0
1130  REAL*8 :: stpmax
1131  LOGICAL :: check
1132 
1133  REAL*8, PARAMETER :: TOLF=1.d-10 , tolmin=1.d-6
1134  REAL*8 :: TOLX
1135 
1136  REAL*8 :: qpj(n_vars)
1137 
1138  REAL*8 :: desc_dir2(n_vars)
1139 
1140  REAL*8 :: desc_dir_temp(n_vars)
1141 
1142  normalize_q = .true.
1143  normalize_f = .false.
1144  opt_search_nl = .true.
1145 
1146  coeff_f(1:n_eqns) = 1.d0
1147 
1148  qj_init = qj
1149 
1150  ! normalize the functions of the nonlinear system
1151  IF ( normalize_f ) THEN
1152 
1153  qj = qj_old - dt * ( matmul(divfluxj+ expl_terms_j,a_tilde) &
1154  - matmul(nhj,a_dirk) )
1155 
1156  CALL eval_f( qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , &
1157  right_term , scal_f )
1158 
1159  IF ( verbose_level .GE. 3 ) THEN
1160 
1161  WRITE(*,*) 'solve_rk_step: non-normalized right_term'
1162  WRITE(*,*) right_term
1163  WRITE(*,*) 'scal_f',scal_f
1164 
1165  END IF
1166 
1167  DO i=1,n_eqns
1168 
1169  IF ( abs(right_term(i)) .GE. 1.d0 ) coeff_f(i) = 1.d0 / right_term(i)
1170 
1171  END DO
1172 
1173  right_term = coeff_f * right_term
1174 
1175  scal_f = 0.5d0 * dot_product( right_term , right_term )
1176 
1177  IF ( verbose_level .GE. 3 ) THEN
1178  WRITE(*,*) 'solve_rk_step: after normalization',scal_f
1179  END IF
1180 
1181  END IF
1182 
1183  !---- normalize the conservative variables ------
1184 
1185  IF ( normalize_q ) THEN
1186 
1187  qj_org = qj
1188 
1189  qj_org = max( abs(qj_org) , 1.d-3 )
1190 
1191  ELSE
1192 
1193  qj_org(1:n_vars) = 1.d0
1194 
1195  END IF
1196 
1197  qj_rel = qj / qj_org
1198 
1199  ! -----------------------------------------------
1200 
1201  newton_raphson_loop:DO nl_iter=1,max_nl_iter
1202 
1203  tolx = epsilon(qj_rel)
1204 
1205  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'solve_rk_step: nl_iter',nl_iter
1206 
1207  CALL eval_f( qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , &
1208  right_term , scal_f )
1209 
1210  IF ( verbose_level .GE. 2 ) THEN
1211 
1212  WRITE(*,*) 'solve_rk_step: right_term',right_term
1213 
1214  END IF
1215 
1216  IF ( verbose_level .GE. 2 ) THEN
1217 
1218  WRITE(*,*) 'before_lnsrch: scal_f',scal_f
1219 
1220  END IF
1221 
1222  ! check the residual of the system
1223 
1224  IF ( maxval( abs( right_term(:) ) ) < tolf ) THEN
1225 
1226  IF ( verbose_level .GE. 3 ) WRITE(*,*) '1: check',check
1227  RETURN
1228 
1229  END IF
1230 
1231  IF ( ( normalize_f ) .AND. ( scal_f < 1.d-6 ) ) THEN
1232 
1233  IF ( verbose_level .GE. 3 ) WRITE(*,*) 'check scal_f',check
1234  RETURN
1235 
1236  END IF
1237 
1238  ! ---- evaluate the descent direction ------------------------------------
1239 
1240  IF ( count( implicit_flag ) .EQ. n_eqns ) THEN
1241 
1242  CALL eval_jacobian( qj_rel , qj_org,coeff_f , left_matrix )
1243 
1244  desc_dir_temp = - right_term
1245 
1246  CALL dgesv(n_eqns,1, left_matrix , n_eqns, pivot, desc_dir_temp , &
1247  n_eqns, ok)
1248 
1249  desc_dir = desc_dir_temp
1250 
1251  ELSE
1252 
1253  CALL eval_jacobian( qj_rel , qj_org,coeff_f , left_matrix )
1254 
1255  left_matrix_small11 = reshape(pack(left_matrix, mask11), &
1256  [n_eqns-n_nh,n_eqns-n_nh])
1257 
1258  left_matrix_small12 = reshape(pack(left_matrix, mask12), &
1259  [n_nh,n_eqns-n_nh])
1260 
1261  left_matrix_small22 = reshape(pack(left_matrix, mask22), &
1262  [n_nh,n_nh])
1263 
1264  left_matrix_small21 = reshape(pack(left_matrix, mask21), &
1265  [n_eqns-n_nh,n_nh])
1266 
1267 
1268  desc_dir_small1 = pack( right_term, .NOT.implicit_flag )
1269  desc_dir_small2 = pack( right_term , implicit_flag )
1270 
1271  DO i=1,n_vars-n_nh
1272 
1273  desc_dir_small1(i) = desc_dir_small1(i) / left_matrix_small11(i,i)
1274 
1275  END DO
1276 
1277  desc_dir_small2 = desc_dir_small2 - &
1278  matmul( desc_dir_small1 , left_matrix_small21 )
1279 
1280  CALL dgesv(n_nh,1, left_matrix_small22 , n_nh , pivot_small2 , &
1281  desc_dir_small2 , n_nh, ok)
1282 
1283  desc_dir = unpack( - desc_dir_small2 , implicit_flag , 0.0d0 ) &
1284  + unpack( - desc_dir_small1 , .NOT.implicit_flag , 0.0d0 )
1285 
1286  END IF
1287 
1288  mod_desc_dir = dsqrt( desc_dir(2)**2 + desc_dir(3)**2 )
1289 
1290  !IF ( qj(2)**2 + qj(3)**2 .GT. 0.D0 ) THEN
1291  !
1292  ! desc_dir(2) = mod_desc_dir * qj(2) / ( qj(2)**2 + qj(3)**2 )
1293  ! desc_dir(3) = mod_desc_dir * qj(3) / ( qj(2)**2 + qj(3)**2 )
1294  !
1295  !END IF
1296 
1297  IF ( verbose_level .GE. 3 ) WRITE(*,*) 'desc_dir',desc_dir
1298 
1299  qj_rel_nr_old = qj_rel
1300  scal_f_old = scal_f
1301 
1302  IF ( ( opt_search_nl ) .AND. ( nl_iter .GT. 1 ) ) THEN
1303  ! Search for the step lambda giving a suffic. decrease in the solution
1304 
1305  stpmax = stpmx * max( sqrt( dot_product(qj_rel,qj_rel) ) , &
1306  dble( SIZE(qj_rel) ) )
1307 
1308  grad_f = matmul( right_term , left_matrix )
1309 
1310  desc_dir2 = desc_dir
1311 
1312  CALL lnsrch( qj_rel_nr_old , qj_org , qj_old , scal_f_old , grad_f , &
1313  desc_dir , coeff_f , qj_rel , scal_f , right_term , stpmax , &
1314  check )
1315 
1316  ELSE
1317 
1318  qj_rel = qj_rel_nr_old + desc_dir
1319 
1320  qj = qj_rel * qj_org
1321 
1322  CALL eval_f( qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , &
1323  right_term , scal_f )
1324 
1325  END IF
1326 
1327  IF ( verbose_level .GE. 2 ) WRITE(*,*) 'after_lnsrch: scal_f',scal_f
1328 
1329  qj = qj_rel * qj_org
1330 
1331  IF ( verbose_level .GE. 3 ) THEN
1332 
1333  WRITE(*,*) 'qj',qj
1334  CALL qc_to_qp( qj , bj , qpj)
1335  WRITE(*,*) 'qp',qpj
1336 
1337  END IF
1338 
1339 
1340  IF ( maxval( abs( right_term(:) ) ) < tolf ) THEN
1341 
1342  IF ( verbose_level .GE. 3 ) WRITE(*,*) '1: check',check
1343  check= .false.
1344  RETURN
1345 
1346  END IF
1347 
1348  IF (check) THEN
1349 
1350  check = ( maxval( abs(grad_f(:)) * max( abs( qj_rel(:) ),1.d0 ) / &
1351  max( scal_f , 0.5d0 * SIZE(qj_rel) ) ) < tolmin )
1352 
1353  IF ( verbose_level .GE. 3 ) WRITE(*,*) '2: check',check
1354  ! RETURN
1355 
1356  END IF
1357 
1358  IF ( maxval( abs( qj_rel(:) - qj_rel_nr_old(:) ) / max( abs( qj_rel(:)) ,&
1359  1.d0 ) ) < tolx ) THEN
1360 
1361  IF ( verbose_level .GE. 3 ) WRITE(*,*) 'check',check
1362  RETURN
1363 
1364  END IF
1365 
1366  END DO newton_raphson_loop
1367 
1368  END SUBROUTINE solve_rk_step
1369 
1370  !******************************************************************************
1372  !
1393  !******************************************************************************
1394 
1395  SUBROUTINE lnsrch( qj_rel_NR_old , qj_org , qj_old , scal_f_old , grad_f , &
1396  desc_dir , coeff_f , qj_rel , scal_f , right_term , stpmax , check )
1397 
1398  IMPLICIT NONE
1399 
1401  REAL*8, DIMENSION(:), INTENT(IN) :: qj_rel_NR_old
1402 
1404  REAL*8, DIMENSION(:), INTENT(IN) :: qj_org
1405 
1407  REAL*8, DIMENSION(:), INTENT(IN) :: qj_old
1408 
1410  REAL*8, DIMENSION(:), INTENT(IN) :: grad_f
1411 
1413  REAL*8, INTENT(IN) :: scal_f_old
1414 
1416  REAL*8, DIMENSION(:), INTENT(INOUT) :: desc_dir
1417 
1418  REAL*8, INTENT(IN) :: stpmax
1419 
1421  REAL*8, DIMENSION(:), INTENT(IN) :: coeff_f
1422 
1424  REAL*8, DIMENSION(:), INTENT(OUT) :: qj_rel
1425 
1427  REAL*8, INTENT(OUT) :: scal_f
1428 
1430  REAL*8, INTENT(OUT) :: right_term(n_eqns)
1431 
1433  LOGICAL, INTENT(OUT) :: check
1434 
1435  REAL*8, PARAMETER :: TOLX=epsilon(qj_rel)
1436 
1437  INTEGER, DIMENSION(1) :: ndum
1438  REAL*8 :: ALF , a,alam,alam2,alamin,b,disc
1439  REAL*8 :: scal_f2
1440  REAL*8 :: desc_dir_abs
1441  REAL*8 :: rhs1 , rhs2 , slope, tmplam
1442 
1443  REAL*8 :: scal_f_min , alam_min
1444 
1445  REAL*8 :: qj(n_vars)
1446 
1447  alf = 1.0d-4
1448 
1449  IF ( size(grad_f) == size(desc_dir) .AND. size(grad_f) == size(qj_rel) &
1450  .AND. size(qj_rel) == size(qj_rel_nr_old) ) THEN
1451 
1452  ndum = size(grad_f)
1453 
1454  ELSE
1455 
1456  WRITE(*,*) 'nrerror: an assert_eq failed with this tag:', 'lnsrch'
1457  stop 'program terminated by assert_eq4'
1458 
1459  END IF
1460 
1461  check = .false.
1462 
1463  desc_dir_abs = sqrt( dot_product(desc_dir,desc_dir) )
1464 
1465  IF ( desc_dir_abs > stpmax ) desc_dir(:) = desc_dir(:) * stpmax/desc_dir_abs
1466 
1467  slope = dot_product(grad_f,desc_dir)
1468 
1469  alamin = tolx / maxval( abs( desc_dir(:))/max( abs(qj_rel_nr_old(:)),1.d0 ) )
1470 
1471  IF ( alamin .EQ. 0.d0) THEN
1472 
1473  qj_rel(:) = qj_rel_nr_old(:)
1474 
1475  RETURN
1476 
1477  END IF
1478 
1479  alam = 1.0d0
1480 
1481  scal_f_min = scal_f_old
1482 
1483  optimal_step_search: DO
1484 
1485  IF ( verbose_level .GE. 4 ) THEN
1486 
1487  WRITE(*,*) 'alam',alam
1488 
1489  END IF
1490 
1491  qj_rel = qj_rel_nr_old + alam * desc_dir
1492 
1493  qj = qj_rel * qj_org
1494 
1495  CALL eval_f( qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , &
1496  right_term , scal_f )
1497 
1498  IF ( verbose_level .GE. 4 ) THEN
1499 
1500  WRITE(*,*) 'lnsrch: effe_old,effe',scal_f_old,scal_f
1501  READ(*,*)
1502 
1503  END IF
1504 
1505  IF ( scal_f .LT. scal_f_min ) THEN
1506 
1507  scal_f_min = scal_f
1508  alam_min = alam
1509 
1510  END IF
1511 
1512  IF ( scal_f .LE. 0.9 * scal_f_old ) THEN
1513  ! sufficient function decrease
1514 
1515  IF ( verbose_level .GE. 4 ) THEN
1516 
1517  WRITE(*,*) 'sufficient function decrease'
1518 
1519  END IF
1520 
1521  EXIT optimal_step_search
1522 
1523  ELSE IF ( alam < alamin ) THEN
1524  ! convergence on Delta_x
1525 
1526  IF ( verbose_level .GE. 4 ) THEN
1527 
1528  WRITE(*,*) ' convergence on Delta_x',alam,alamin
1529 
1530  END IF
1531 
1532  qj_rel(:) = qj_rel_nr_old(:)
1533  scal_f = scal_f_old
1534  check = .true.
1535 
1536  EXIT optimal_step_search
1537 
1538  ! ELSE IF ( scal_f .LE. scal_f_old + ALF * alam * slope ) THEN
1539  ELSE
1540 
1541  IF ( alam .EQ. 1.d0 ) THEN
1542 
1543  tmplam = - slope / ( 2.0d0 * ( scal_f - scal_f_old - slope ) )
1544 
1545  ELSE
1546 
1547  rhs1 = scal_f - scal_f_old - alam*slope
1548  rhs2 = scal_f2 - scal_f_old - alam2*slope
1549 
1550  a = ( rhs1/alam**2.d0 - rhs2/alam2**2.d0 ) / ( alam - alam2 )
1551  b = ( -alam2*rhs1/alam**2 + alam*rhs2/alam2**2 ) / ( alam - alam2 )
1552 
1553  IF ( a .EQ. 0.d0 ) THEN
1554 
1555  tmplam = - slope / ( 2.0d0 * b )
1556 
1557  ELSE
1558 
1559  disc = b*b - 3.0d0*a*slope
1560 
1561  IF ( disc .LT. 0.d0 ) THEN
1562 
1563  tmplam = 0.5d0 * alam
1564 
1565  ELSE IF ( b .LE. 0.d0 ) THEN
1566 
1567  tmplam = ( - b + sqrt(disc) ) / ( 3.d0 * a )
1568 
1569  ELSE
1570 
1571  tmplam = - slope / ( b + sqrt(disc) )
1572 
1573  ENDIF
1574 
1575  END IF
1576 
1577  IF ( tmplam .GT. 0.5d0*alam ) tmplam = 0.5d0 * alam
1578 
1579  END IF
1580 
1581  END IF
1582 
1583  alam2 = alam
1584  scal_f2 = scal_f
1585  alam = max( tmplam , 0.5d0*alam )
1586 
1587  END DO optimal_step_search
1588 
1589  END SUBROUTINE lnsrch
1590 
1591  !******************************************************************************
1593  !
1609  !******************************************************************************
1610 
1611  SUBROUTINE eval_f( qj , qj_old , a_tilde , a_dirk , a_diag , coeff_f , f_nl , &
1612  scal_f )
1613 
1614  USE constitutive_2d, ONLY : eval_nonhyperbolic_terms
1615 
1616  IMPLICIT NONE
1617 
1618  REAL*8, INTENT(IN) :: qj(n_vars)
1619  REAL*8, INTENT(IN) :: qj_old(n_vars)
1620  REAL*8, INTENT(IN) :: a_tilde(n_rk)
1621  REAL*8, INTENT(IN) :: a_dirk(n_rk)
1622  REAL*8, INTENT(IN) :: a_diag
1623  REAL*8, INTENT(IN) :: coeff_f(n_eqns)
1624  REAL*8, INTENT(OUT) :: f_nl(n_eqns)
1625  REAL*8, INTENT(OUT) :: scal_f
1626 
1627  REAL*8 :: nh_term_impl(n_eqns)
1628  REAL*8 :: Rj(n_eqns)
1629 
1630  CALL eval_nonhyperbolic_terms( r_qj = qj , r_nh_term_impl=nh_term_impl )
1631 
1632  rj = ( matmul( divfluxj(1:n_eqns,1:i_rk-1)+expl_terms_j(1:n_eqns,1:i_rk-1), &
1633  a_tilde(1:i_rk-1) ) - matmul( nhj(1:n_eqns,1:i_rk-1) &
1634  + si_nhj(1:n_eqns,1:i_rk-1) , a_dirk(1:i_rk-1) ) ) &
1635  - a_diag * ( nh_term_impl + si_nhj(1:n_eqns,i_rk) )
1636 
1637  f_nl = qj - qj_old + dt * rj
1638 
1639  f_nl = coeff_f * f_nl
1640 
1641  scal_f = 0.5d0 * dot_product( f_nl , f_nl )
1642 
1643  !WRITE(*,*) 'i_RK',i_RK
1644  !WRITE(*,*) 'eval_f: qj(2) = ',qj(2)
1645  !WRITE(*,*) 'nh_term_impl(2) = ',nh_term_impl(2)
1646  !WRITE(*,*) ' SI_NHj(1:n_eqns,i_RK) = ', SI_NHj(1:n_eqns,i_RK)
1647  !WRITE(*,*) 'eval_f: new term = ',qj(2) + dt * ( - a_diag * nh_term_impl(2) )
1648  !WRITE(*,*) 'f_nl(2)=',f_nl(2)
1649  !READ(*,*)
1650 
1651  END SUBROUTINE eval_f
1652 
1653  !******************************************************************************
1655  !
1658  !
1667  !
1671  !******************************************************************************
1672 
1673  SUBROUTINE eval_jacobian( qj_rel , qj_org , coeff_f, left_matrix)
1674 
1675  USE constitutive_2d, ONLY : eval_nonhyperbolic_terms
1676 
1677  IMPLICIT NONE
1678 
1679  REAL*8, INTENT(IN) :: qj_rel(n_vars)
1680  REAL*8, INTENT(IN) :: qj_org(n_vars)
1681  REAL*8, INTENT(IN) :: coeff_f(n_eqns)
1682 
1683  REAL*8, INTENT(OUT) :: left_matrix(n_eqns,n_vars)
1684 
1685  REAL*8 :: Jacob_relax(n_eqns,n_vars)
1686  COMPLEX*16 :: nh_terms_cmplx_impl(n_eqns)
1687  COMPLEX*16 :: qj_cmplx(n_vars) , qj_rel_cmplx(n_vars)
1688 
1689  REAL*8 :: h
1690 
1691  INTEGER :: i
1692 
1693  h = n_vars * epsilon(1.d0)
1694 
1695  ! initialize the matrix of the linearized system and the Jacobian
1696 
1697  left_matrix(1:n_eqns,1:n_vars) = 0.d0
1698  jacob_relax(1:n_eqns,1:n_vars) = 0.d0
1699 
1700  ! evaluate the jacobian of the non-hyperbolic terms
1701 
1702  DO i=1,n_vars
1703 
1704  left_matrix(i,i) = coeff_f(i) * qj_org(i)
1705 
1706  IF ( implicit_flag(i) ) THEN
1707 
1708  qj_rel_cmplx(1:n_vars) = qj_rel(1:n_vars)
1709  qj_rel_cmplx(i) = dcmplx(qj_rel(i), h)
1710 
1711  qj_cmplx = qj_rel_cmplx * qj_org
1712 
1713  CALL eval_nonhyperbolic_terms( c_qj = qj_cmplx , &
1714  c_nh_term_impl = nh_terms_cmplx_impl )
1715 
1716  jacob_relax(1:n_eqns,i) = coeff_f(i) * &
1717  aimag(nh_terms_cmplx_impl) / h
1718 
1719  left_matrix(1:n_eqns,i) = left_matrix(1:n_eqns,i) - dt * a_diag &
1720  * jacob_relax(1:n_eqns,i)
1721 
1722  END IF
1723 
1724  END DO
1725 
1726  END SUBROUTINE eval_jacobian
1727 
1728 
1729  !******************************************************************************
1731  !
1735  !
1738  !
1742  !******************************************************************************
1743 
1744  SUBROUTINE update_erosion_deposition_ver(dt)
1745 
1746  USE constitutive_2d, ONLY : eval_erosion_dep_term
1747  USE constitutive_2d, ONLY : eval_topo_term
1748 
1749  USE constitutive_2d, ONLY : qc_to_qp
1750 
1751  USE constitutive_2d, ONLY : erosion_coeff , settling_vel
1752  USE geometry_2d, ONLY : b_nw , b_ne , b_sw , b_se
1753 
1754  IMPLICIT NONE
1755 
1756  REAL*8, INTENT(IN) :: dt
1757 
1758  REAL*8 :: erosion_termLT , erosion_termRT
1759  REAL*8 :: erosion_termLB , erosion_termRB
1760 
1761  REAL*8 :: deposition_termLT , deposition_termRT
1762  REAL*8 :: deposition_termLB , deposition_termRB
1763 
1764  REAL*8 :: vertex_erosion_terms(comp_interfaces_x,comp_interfaces_y)
1765  REAL*8 :: vertex_deposition_terms(comp_interfaces_x,comp_interfaces_y)
1766 
1767 
1768  REAL*8 :: avg_erosion_term
1769  REAL*8 :: avg_deposition_term
1770  REAL*8 :: eqns_term(n_eqns)
1771  REAL*8 :: topo_term
1772 
1773  REAL*8 :: B_avg
1774 
1775  INTEGER :: j ,k
1776 
1777  IF ( ( erosion_coeff .EQ. 0.d0 ) .AND. ( settling_vel .EQ. 0.d0 ) ) RETURN
1778 
1779  ! Bi-linear reconstruction of the solution at the volume interfaces
1780  CALL reconstruction
1781 
1782  ! Loop to extend the bilinear reconstruction at the 4 corners of each volume
1783  DO k = 1,comp_cells_y
1784 
1785  DO j = 1,comp_cells_x
1786 
1787  ! Bilinear reconstruction of the conservative variables at the corners
1788  q_cellnw(:,j,k) = q(:,j,k) - 0.5d0 * ( q_interfacel(:,j+1,k) - &
1789  q_interfacer(:,j,k) ) + 0.5d0 * ( q_interfaceb(:,j,k+1) - &
1790  q_interfacet(:,j,k) )
1791 
1792  q_cellne(:,j,k) = q(:,j,k) + 0.5d0 * ( q_interfacel(:,j+1,k) - &
1793  q_interfacer(:,j,k) ) + 0.5d0 * ( q_interfaceb(:,j,k+1) - &
1794  q_interfacet(:,j,k) )
1795 
1796  q_cellsw(:,j,k) = q(:,j,k) - 0.5d0 * ( q_interfacel(:,j+1,k) - &
1797  q_interfacer(:,j,k) ) - 0.5d0 * ( q_interfaceb(:,j,k+1) - &
1798  q_interfacet(:,j,k) )
1799 
1800  q_cellse(:,j,k) = q(:,j,k) + 0.5d0 * ( q_interfacel(:,j+1,k) - &
1801  q_interfacer(:,j,k) ) - 0.5d0 * ( q_interfaceb(:,j,k+1) - &
1802  q_interfacet(:,j,k) )
1803 
1804  ! Bilinear reconstruction of the topography at the corners
1805  b_nw(j,k) = b_cent(j,k) - 0.5d0 * ( b_stag_x(j+1,k) - b_stag_x(j,k) ) &
1806  + 0.5d0 * ( b_stag_y(j,k+1) - b_stag_y(j,k) )
1807 
1808  b_ne(j,k) = b_cent(j,k) + 0.5d0 * ( b_stag_x(j+1,k) - b_stag_x(j,k) ) &
1809  + 0.5d0 * ( b_stag_y(j,k+1) - b_stag_y(j,k) )
1810 
1811  b_sw(j,k) = b_cent(j,k) - 0.5d0 * ( b_stag_x(j+1,k) - b_stag_x(j,k) ) &
1812  - 0.5d0 * ( b_stag_y(j,k+1) - b_stag_y(j,k) )
1813 
1814  b_se(j,k) = b_cent(j,k) + 0.5d0 * ( b_stag_x(j+1,k) - b_stag_x(j,k) ) &
1815  - 0.5d0 * ( b_stag_y(j,k+1) - b_stag_y(j,k) )
1816 
1817 
1818  ! Correct the slope of the reconstruction if h<0 at NW corner
1819  IF ( q_cellnw(1,j,k) .LT. 0.d0 ) THEN
1820 
1821  q_cellnw(1,j,k) = 0.d0
1822  q_cellse(1,j,k) = max(0.d0,2.d0 * q(1,j,k))
1823 
1824  q_cellnw(4,j,k) = 0.d0
1825  q_cellse(4,j,k) = max(0.d0,2.d0 * q(4,j,k))
1826 
1827  ENDIF
1828 
1829  ! Correct the slope of the reconstruction if h<0 at SE corner
1830  IF ( q_cellse(1,j,k) .LT. 0.d0 ) THEN
1831 
1832  q_cellse(1,j,k) = 0.d0
1833  q_cellnw(1,j,k) = max(0.d0,2.d0 * q(1,j,k))
1834 
1835  q_cellse(4,j,k) = 0.d0
1836  q_cellnw(4,j,k) = max(0.d0,2.d0 * q(4,j,k))
1837 
1838  ENDIF
1839 
1840  ! Correct the slope of the reconstruction if h<0 at NE corner
1841  IF ( q_cellne(1,j,k) .LT. 0.d0 ) THEN
1842 
1843  q_cellne(1,j,k) = 0.d0
1844  q_cellsw(1,j,k) = max(0.d0,2.d0 * q(1,j,k))
1845 
1846  q_cellne(4,j,k) = 0.d0
1847  q_cellsw(4,j,k) = max(0.d0,2.d0 * q(4,j,k))
1848 
1849  ENDIF
1850 
1851  ! Correct the slope of the reconstruction if h<0 at SW corner
1852  IF ( q_cellsw(1,j,k) .LT. 0.d0 ) THEN
1853 
1854  q_cellsw(1,j,k) = 0.d0
1855  q_cellne(1,j,k) = max(0.d0,2.d0 * q(1,j,k))
1856 
1857  q_cellsw(4,j,k) = 0.d0
1858  q_cellne(4,j,k) = max(0.d0,2.d0 * q(4,j,k))
1859 
1860  ENDIF
1861 
1862 
1863  !-----
1864 
1865  ! Correct the slope of the reconstruction if h*alpha<0 at NW corner
1866  IF ( q_cellnw(4,j,k) .LT.0.d0 ) THEN
1867 
1868  q_cellnw(4,j,k) = 0.d0
1869  q_cellse(4,j,k) = max(0.d0,2.d0 * q(4,j,k))
1870 
1871  ! Correct the slope of the reconstruction if alpha>1 at NW corner
1872  ELSEIF ( q_cellnw(4,j,k) .GT. q_cellnw(1,j,k) ) THEN
1873 
1874  q_cellnw(4,j,k) = q_cellnw(1,j,k)
1875  q_cellse(4,j,k) = 2.d0 * q(4,j,k) - q_cellnw(4,j,k)
1876 
1877  ENDIF
1878 
1879  ! Correct the slope of the reconstruction if h*alpha<0 at SE corner
1880  IF ( q_cellse(4,j,k) .LT. 0.d0 ) THEN
1881 
1882  q_cellse(4,j,k) = 0.d0
1883  q_cellnw(4,j,k) = max(0.d0,2.d0 * q(4,j,k))
1884 
1885  ! Correct the slope of the reconstruction if alpha>1 at SE corner
1886  ELSEIF ( q_cellse(4,j,k) .GT. q_cellse(1,j,k) ) THEN
1887 
1888  q_cellse(4,j,k) = q_cellse(1,j,k)
1889  q_cellnw(4,j,k) = 2.d0 * q(4,j,k) - q_cellse(4,j,k)
1890 
1891  ENDIF
1892 
1893  ! Correct the slope of the reconstruction if h*alpha<0 at NE corner
1894  IF ( q_cellne(4,j,k) .LT.0.d0 ) THEN
1895 
1896  q_cellne(4,j,k) = 0.d0
1897  q_cellsw(4,j,k) = max(0.d0,2.d0 * q(4,j,k))
1898 
1899  ! Correct the slope of the reconstruction if alpha>1 at NE corner
1900  ELSEIF ( q_cellne(4,j,k) .GT. q_cellne(1,j,k) ) THEN
1901 
1902  q_cellne(4,j,k) = q_cellne(1,j,k)
1903  q_cellsw(4,j,k) = 2.d0 * q(4,j,k) - q_cellne(4,j,k)
1904 
1905  ENDIF
1906 
1907  ! Correct the slope of the reconstruction if h*alpha<0 at SW corner
1908  IF ( q_cellsw(4,j,k) .LT. 0.d0 ) THEN
1909 
1910  q_cellsw(4,j,k) = 0.d0
1911  q_cellne(4,j,k) = max(0.d0,2.d0 * q(4,j,k))
1912 
1913  ! Correct the slope of the reconstruction if alpha>1 at SW corner
1914  ELSEIF ( q_cellsw(4,j,k) .GT. q_cellsw(1,j,k) ) THEN
1915 
1916  q_cellsw(4,j,k) = q_cellsw(1,j,k)
1917  q_cellne(4,j,k) = 2.d0 * q(4,j,k) - q_cellsw(4,j,k)
1918 
1919  ENDIF
1920 
1921  ! Final checks and correctionsfor negative values at NW corner
1922  IF ( q_cellnw(4,j,k) .LT. 0.d0 ) THEN
1923 
1924  IF ( q_cellnw(4,j,k) .GT. -1.d-10 ) THEN
1925 
1926  q_cellnw(4,j,k) = 0.d0
1927 
1928  ELSE
1929 
1930  WRITE(*,*) 'j,k',j,k
1931  WRITE(*,*) 'q_cellNW(4,j,k)', q_cellnw(4,j,k)
1932  WRITE(*,*) 'q_cellSE(4,j,k)', q_cellse(4,j,k)
1933  WRITE(*,*) 'q(4,j,k)',q(4,j,k)
1934  READ(*,*)
1935 
1936  END IF
1937 
1938  END IF
1939 
1940  ! Final checks and correctionsfor negative values at NE corner
1941  IF ( q_cellne(4,j,k) .LT. 0.d0 ) THEN
1942 
1943  IF ( q_cellne(4,j,k) .GT. -1.d-10 ) THEN
1944 
1945  q_cellne(4,j,k) = 0.d0
1946 
1947  ELSE
1948 
1949  WRITE(*,*) 'j,k',j,k
1950  WRITE(*,*) ' q_cellNE(4,j,k)', q_cellne(4,j,k)
1951  READ(*,*)
1952 
1953  END IF
1954 
1955  END IF
1956 
1957  ! Final checks and correctionsfor negative values at SW corner
1958  IF ( q_cellsw(4,j,k) .LT. 0.d0 ) THEN
1959 
1960  IF ( q_cellsw(4,j,k) .GT. -1.d-10 ) THEN
1961 
1962  q_cellsw(4,j,k) = 0.d0
1963 
1964  ELSE
1965 
1966  WRITE(*,*) 'j,k',j,k
1967  WRITE(*,*) ' q_cellSW(4,j,k)', q_cellsw(4,j,k)
1968  WRITE(*,*) q_interfacer(4,j,k)
1969  WRITE(*,*) q_interfacel(4,j+1,k)
1970  WRITE(*,*) q_interfacet(4,j,k)
1971  WRITE(*,*) q_interfaceb(4,j,k+1)
1972  WRITE(*,*) q(4,j,k)
1973  READ(*,*)
1974 
1975  END IF
1976 
1977  END IF
1978 
1979  ! Final checks and correctionsfor negative values at SE corner
1980  IF ( q_cellse(4,j,k) .LT. 0.d0 ) THEN
1981 
1982  IF ( q_cellse(4,j,k) .GT. -1.d-10 ) THEN
1983 
1984  q_cellse(4,j,k) = 0.d0
1985 
1986  ELSE
1987 
1988  WRITE(*,*) 'j,k',j,k
1989  WRITE(*,*) ' q_cellSE(4,j,k)', q_cellse(4,j,k)
1990  READ(*,*)
1991 
1992  END IF
1993 
1994  END IF
1995 
1996  END DO
1997 
1998  END DO
1999 
2000  ! Initialization of erosion and deposition at vertices of the cells
2001  vertex_erosion_terms(1:comp_interfaces_x,1:comp_interfaces_y) = 0.d0
2002  vertex_deposition_terms(1:comp_interfaces_x,1:comp_interfaces_y) = 0.d0
2003 
2004  ! ---------------------------------------------------------------------------
2005  ! Compute the erosion and deposition terms at the south-boundary vertexes
2006  ! of the domain (k=1)
2007 
2008  k = 1
2009  j = 1
2010 
2011  CALL eval_erosion_dep_term( q_cellsw(:,j,k) , erosion_termrt , &
2012  deposition_termrt )
2013 
2014  deposition_termrt = min( deposition_termrt , q_cellsw(4,j,k) / dt )
2015 
2016  vertex_erosion_terms(j,k) = erosion_termrt
2017 
2018  vertex_deposition_terms(j,k) = deposition_termrt
2019 
2020  south_loop:DO j = 2,comp_interfaces_x-1
2021 
2022  ! ------- Erosion and deposition at the SW side of vertex ----------------
2023  CALL eval_erosion_dep_term( q_cellsw(:,j,k) , erosion_termrt , &
2024  deposition_termrt )
2025 
2026  deposition_termrt = min( deposition_termrt , q_cellsw(4,j,k) / dt )
2027 
2028  ! ------- Erosion and deposition at the SE side of vertex ----------------
2029  CALL eval_erosion_dep_term( q_cellse(:,j-1,k) , erosion_termlt , &
2030  deposition_termlt )
2031 
2032  deposition_termlt = min( deposition_termlt , q_cellse(4,j-1,k) / dt )
2033 
2034  ! ------------------- Erosion at the vertex (j,k) ------------------------
2035  vertex_erosion_terms(j,k) = min( erosion_termrt , erosion_termlt )
2036 
2037  ! ----------------- Deposition at the vertex (j,k) -----------------------
2038  vertex_deposition_terms(j,k) = min( deposition_termrt, deposition_termlt )
2039 
2040  END DO south_loop
2041 
2042  j = comp_interfaces_x
2043 
2044  ! ------- Erosion and deposition at the SE side of vertex -------------------
2045  CALL eval_erosion_dep_term( q_cellse(:,j-1,k) , erosion_termlt , &
2046  deposition_termlt )
2047 
2048  deposition_termlt = min( deposition_termlt , q_cellse(4,j-1,k) / dt )
2049 
2050  vertex_erosion_terms(j,k) = erosion_termlt
2051 
2052  vertex_deposition_terms(j,k) = deposition_termlt
2053 
2054  ! End loop for erosion and deposition terms at the south-boundary vertexes
2055  ! of the domain (k=1)
2056  ! ---------------------------------------------------------------------------
2057 
2058 
2059  ! ---------------------------------------------------------------------------
2060  ! Compute the erosion and deposition terms at the west-boundary vertexes
2061  ! of the domain, corners excluded (j=1, k=2,comp_interfaces_y-1)
2062 
2063  j = 1
2064 
2065  west_loop:DO k = 2,comp_interfaces_y-1
2066 
2067  ! ------- Erosion and deposition at the SW side of vertex ----------------
2068  CALL eval_erosion_dep_term( q_cellsw(:,j,k) , erosion_termrt , &
2069  deposition_termrt )
2070 
2071  deposition_termrt = min( deposition_termrt , q_cellsw(4,j,k) / dt )
2072 
2073  ! ------- Erosion and deposition at the NW side of vertex ----------------
2074  CALL eval_erosion_dep_term( q_cellnw(:,j,k-1) , erosion_termrb , &
2075  deposition_termrb )
2076 
2077  deposition_termrb = min( deposition_termrb , q_cellnw(4,j,k-1) / dt )
2078 
2079  ! ------------------- Erosion at the (j,k) vertex ------------------------
2080  vertex_erosion_terms(j,k) = min( erosion_termrt , erosion_termrb )
2081 
2082  ! ----------------- Deposition at the (j,k) vertex -----------------------
2083  vertex_deposition_terms(j,k) = min( deposition_termrt, deposition_termrb )
2084 
2085  END DO west_loop
2086 
2087  ! End loop for erosion and deposition terms at the west-boundary vertexes
2088  ! of the domain, corners excluded (j=1, k=2,comp_interfaces_y-1)
2089  ! ---------------------------------------------------------------------------
2090 
2091  ! ---------------------------------------------------------------------------
2092  ! Compute the erosion and deposition terms at the east-boundary vertexes of
2093  ! the domain, corners excluded (j=comp_interfaces_x, k=2,comp_interfaces_y-1)
2094 
2095  j = comp_interfaces_x
2096 
2097  east_loop:DO k = 2,comp_interfaces_y-1
2098 
2099  ! ------- Erosion and deposition at the NE side of vertex ----------------
2100  CALL eval_erosion_dep_term( q_cellne(:,j-1,k-1) , erosion_termlb , &
2101  deposition_termlb )
2102 
2103  deposition_termlb = min( deposition_termlb , q_cellne(4,j-1,k-1) / dt )
2104 
2105  ! ------- Erosion and deposition at the SE side of vertex ----------------
2106  CALL eval_erosion_dep_term( q_cellse(:,j-1,k) , erosion_termlt , &
2107  deposition_termlt )
2108 
2109  deposition_termlt = min( deposition_termlt , q_cellse(4,j-1,k) / dt )
2110 
2111  ! ------------------- Erosion at the (j,k) vertex ------------------------
2112  vertex_erosion_terms(j,k) = min( erosion_termlt , erosion_termlb )
2113 
2114  ! ----------------- Deposition at the (j,k) vertex -----------------------
2115  vertex_deposition_terms(j,k) = min( deposition_termlt , deposition_termlb)
2116 
2117  END DO east_loop
2118 
2119  ! End loop for erosion and deposition terms at the east-boundary vertexes of
2120  ! the domain, corners excluded (j=comp_interfaces_x, k=2,comp_interfaces_y-1)
2121  ! ---------------------------------------------------------------------------
2122 
2123 
2124  ! ---------------------------------------------------------------------------
2125  ! Compute the erosion and deposition terms at the north-boundary vertexes
2126  ! of the domain (k=comp_interfaces_y , j = 1,comp_interfaces_x)
2127 
2128 
2129  ! End loop for erosion and deposition terms at the north-boundary vertexes
2130  ! of the domain (k=comp_interfaces_y , j = 1,comp_interfaces_x)
2131  ! ---------------------------------------------------------------------------
2132 
2133  k = comp_interfaces_y
2134  j = 1
2135 
2136  ! ------- Erosion and deposition at the NW side of vertex -------------------
2137  CALL eval_erosion_dep_term( q_cellnw(:,j,k-1) , erosion_termrb , &
2138  deposition_termrb )
2139 
2140  IF ( erosion_termrb .LT. 0.d0 ) THEN
2141 
2142  WRITE(*,*) 'j,k',j,k
2143  WRITE(*,*) 'erosion_termRB',erosion_termrb
2144  WRITE(*,*) 'q_cellNW(:,j,k-1)',q_cellnw(:,j,k-1)
2145  READ(*,*)
2146 
2147  END IF
2148 
2149  deposition_termrb = min( deposition_termrb , q_cellnw(4,j,k-1) / dt )
2150 
2151  vertex_erosion_terms(j,k) = erosion_termrb
2152  vertex_deposition_terms(j,k) = deposition_termrb
2153 
2154  DO j = 2,comp_interfaces_x-1
2155 
2156  ! ------- Erosion and deposition at the NW side of vertex ----------------
2157  CALL eval_erosion_dep_term( q_cellnw(:,j,k-1) , erosion_termrb , &
2158  deposition_termrb )
2159 
2160  deposition_termrb = min( deposition_termrb , q_cellnw(4,j,k-1) / dt )
2161 
2162  ! ------- Erosion and deposition at the NE side of vertex ----------------
2163  CALL eval_erosion_dep_term( q_cellne(:,j-1,k-1) , erosion_termlb , &
2164  deposition_termlb )
2165 
2166 
2167  deposition_termlb = min( deposition_termlb , q_cellne(4,j-1,k-1) / dt )
2168 
2169  ! ------------------- Erosion at the (j,k) vertex ------------------------
2170  vertex_erosion_terms(j,k) = min( erosion_termrb , erosion_termlb )
2171 
2172  ! ----------------- Deposition at the (j,k) vertex -----------------------
2173  vertex_deposition_terms(j,k) = min( deposition_termrb , deposition_termlb)
2174 
2175  END DO
2176 
2177  j = comp_interfaces_x
2178 
2179  ! ------- Erosion and deposition at the NE side of vertex -------------------
2180  CALL eval_erosion_dep_term( q_cellne(:,j-1,k-1) , erosion_termlb , &
2181  deposition_termlb )
2182 
2183  deposition_termlb = min( deposition_termlb , q_cellne(4,j-1,k-1) / dt )
2184 
2185  vertex_erosion_terms(j,k) = erosion_termlb
2186  vertex_deposition_terms(j,k) = deposition_termlb
2187 
2188  ! ---------------------------------------------------------------------------
2189  ! Compute the erosion and deposition terms at all internal vertexes
2190  ! of the domain (j=2,comp_interfaces_x-1, k=2,comp_interfaces_y-1)
2191 
2192  internal_vertexes_y:DO k = 2,comp_interfaces_y-1
2193 
2194  internal_vertexes_x:DO j = 2,comp_interfaces_x-1
2195 
2196  ! ------- Erosion and deposition at the SW side of vertex -------------
2197  CALL eval_erosion_dep_term( q_cellsw(:,j,k) , erosion_termrt , &
2198  deposition_termrt )
2199 
2200  IF ( erosion_termrt .LT. 0.d0 ) THEN
2201 
2202  WRITE(*,*) 'j,k',j,k
2203  WRITE(*,*) 'erosion_termRT',erosion_termrt
2204  WRITE(*,*) 'q_cellSW(:,j,k)',q_cellsw(:,j,k)
2205  READ(*,*)
2206 
2207  END IF
2208 
2209  deposition_termrt = min( deposition_termrt , q_cellsw(4,j,k) / dt )
2210 
2211  ! ------- Erosion and deposition at the SE side of vertex -------------
2212  CALL eval_erosion_dep_term( q_cellse(:,j-1,k) , erosion_termlt , &
2213  deposition_termlt )
2214 
2215  IF ( erosion_termlt .LT. 0.d0 ) THEN
2216 
2217  WRITE(*,*) 'j,k',j,k
2218  WRITE(*,*) 'erosion_termLT',erosion_termlt
2219  WRITE(*,*) 'q_cellSE(:,j-1,k)',q_cellse(:,j-1,k)
2220  READ(*,*)
2221 
2222  END IF
2223 
2224  deposition_termlt = min( deposition_termlt , q_cellse(4,j-1,k) / dt )
2225 
2226  ! ------- Erosion and deposition at the NW side of vertex -------------
2227  CALL eval_erosion_dep_term( q_cellnw(:,j,k-1) , erosion_termrb , &
2228  deposition_termrb )
2229 
2230  IF ( erosion_termrb .LT. 0.d0 ) THEN
2231 
2232  WRITE(*,*) 'j,k',j,k
2233  WRITE(*,*) 'erosion_termRB',erosion_termrb
2234  WRITE(*,*) 'q_cellNW(:,j,k-1)',q_cellnw(:,j,k-1)
2235  READ(*,*)
2236 
2237  END IF
2238 
2239  deposition_termrb = min( deposition_termrb , q_cellnw(4,j,k-1) / dt )
2240 
2241  ! ------- Erosion and deposition at the NE side of vertex -------------
2242  CALL eval_erosion_dep_term( q_cellne(:,j-1,k-1) , erosion_termlb , &
2243  deposition_termlb )
2244 
2245  IF ( erosion_termlb .LT. 0.d0 ) THEN
2246 
2247  WRITE(*,*) 'j,k',j,k
2248  WRITE(*,*) 'erosion_termLB',erosion_termlb
2249  WRITE(*,*) 'q_cellNE(:,j-1,k-1)',q_cellne(:,j-1,k-1)
2250  READ(*,*)
2251 
2252  END IF
2253 
2254  deposition_termlb = min( deposition_termlb , q_cellne(4,j-1,k-1) / dt )
2255 
2256  ! ------------------- Erosion at the (j,k) vertex ---------------------
2257  vertex_erosion_terms(j,k) = min( erosion_termrt , erosion_termlt , &
2258  erosion_termrb , erosion_termlb )
2259 
2260  IF (isnan(vertex_erosion_terms(j,k))) THEN
2261 
2262  WRITE(*,*) 'j,k',j,k
2263  WRITE(*,*) erosion_termrt , erosion_termlt , erosion_termrb , &
2264  erosion_termlb
2265  WRITE(*,*) b_se(j,k) , q_cellse(:,j-1,k)
2266  stop
2267 
2268  END IF
2269 
2270  ! ----------------- Deposition at the (j,k) vertex --------------------
2271  vertex_deposition_terms(j,k) = min( deposition_termrt , &
2272  deposition_termlt , deposition_termrb , deposition_termlb )
2273 
2274  END DO internal_vertexes_x
2275 
2276  END DO internal_vertexes_y
2277 
2278  ! ---------------------------------------------------------------------------
2279  ! End loop for erosion and deposition terms at all internal vertexes
2280  ! of the domain (j=2,comp_interfaces_x-1, k=2,comp_interfaces_y-1)
2281 
2282 
2283  ! --------------- Update topography at the all vertexes ---------------------
2284  vertexes_y:DO k = 1,comp_interfaces_y
2285 
2286  vertexes_x:DO j = 1,comp_interfaces_x
2287 
2288  b_ver(j,k) = b_ver(j,k) + dt * ( - vertex_erosion_terms(j,k) + &
2289  vertex_deposition_terms(j,k) )
2290 
2291  END DO vertexes_x
2292 
2293  END DO vertexes_y
2294 
2295  ! ------------ Compute erosion and deposition at the cell faces -------------
2296  DO k = 1,comp_interfaces_y
2297 
2298  DO j = 1,comp_cells_x
2299 
2300  avg_erosion_term = 0.5d0 * ( vertex_erosion_terms(j,k) &
2301  + vertex_erosion_terms(j+1,k) )
2302 
2303  avg_deposition_term = 0.5d0 * ( vertex_deposition_terms(j,k) &
2304  + vertex_deposition_terms(j+1,k) )
2305 
2306  b_stag_y(j,k) = b_stag_y(j,k) + dt * ( avg_deposition_term - &
2307  avg_erosion_term )
2308 
2309  END DO
2310 
2311  END DO
2312 
2313  DO k = 1,comp_cells_y
2314 
2315  DO j = 1,comp_interfaces_x
2316 
2317  avg_erosion_term = 0.5d0 * ( vertex_erosion_terms(j,k) &
2318  + vertex_erosion_terms(j,k+1) )
2319 
2320  avg_deposition_term = 0.5d0 * ( vertex_deposition_terms(j,k) &
2321  + vertex_deposition_terms(j,k+1) )
2322 
2323  b_stag_x(j,k) = b_stag_x(j,k) + dt * ( avg_deposition_term - &
2324  avg_erosion_term )
2325 
2326  END DO
2327 
2328  END DO
2329 
2330 
2331  ! For each cell, the erosion and deposition term is obtained with an
2332  ! average of the values computed at the four vertexes of the cell
2333  DO k = 1,comp_cells_y
2334 
2335  DO j = 1,comp_cells_x
2336 
2337  avg_erosion_term = 0.25d0 * ( vertex_erosion_terms(j,k) &
2338  + vertex_erosion_terms(j+1,k) &
2339  + vertex_erosion_terms(j,k+1) &
2340  + vertex_erosion_terms(j+1,k+1) )
2341 
2342  avg_deposition_term = 0.25d0 * ( vertex_deposition_terms(j,k) &
2343  + vertex_deposition_terms(j+1,k) &
2344  + vertex_deposition_terms(j,k+1) &
2345  + vertex_deposition_terms(j+1,k+1) )
2346 
2347  ! Compute the source terms for the equations
2348  CALL eval_topo_term( q(1:n_vars,j,k) , avg_deposition_term , &
2349  avg_erosion_term , eqns_term , topo_term )
2350 
2351  IF ( verbose_level .GE. 2 ) THEN
2352 
2353  WRITE(*,*) 'before update erosion/deposition: j,k,q(:,j,k),B(j,k)',&
2354  j,k,q(:,j,k),b_cent(j,k)
2355 
2356  END IF
2357 
2358  ! Update the topography with erosion/deposition terms
2359  b_cent(j,k) = b_cent(j,k) + dt * topo_term
2360 
2361  ! First check for bi-linearity of recontructed B
2362  IF ( dabs(b_cent(j,k) - 0.5 * ( b_stag_x(j,k) + b_stag_x(j+1,k) ) ) &
2363  .GT. 1.d-5) THEN
2364 
2365  WRITE(*,*) 'after update erosion/deposition'
2366  WRITE(*,*) 'j,k',j,k
2367  WRITE(*,*) 'B val x'
2368  WRITE(*,*) b_cent(j,k) ,0.5*(b_stag_x(j,k) + b_stag_x(j+1,k)) , &
2369  b_stag_x(j,k) , b_stag_x(j+1,k)
2370  READ(*,*)
2371 
2372  END IF
2373 
2374  ! Second check for bi-linearity of recontructed B
2375  IF ( dabs(b_cent(j,k) - 0.5 * ( b_stag_y(j,k) + b_stag_y(j,k+1) ) ) &
2376  .GT. 1.d-5) THEN
2377 
2378  WRITE(*,*) 'after update erosion/deposition'
2379  WRITE(*,*) 'j,k',j,k
2380  WRITE(*,*) 'B val y'
2381  WRITE(*,*) b_cent(j,k) ,0.5*(b_stag_y(j,k) + b_stag_y(j,k+1)) , &
2382  b_stag_y(j,k) , b_stag_y(j,k+1)
2383  READ(*,*)
2384 
2385  END IF
2386 
2387  ! Third check for bi-linearity of recontructed B
2388  b_avg = 0.25d0 * ( b_stag_x(j,k) + b_stag_x(j+1,k ) + b_stag_y(j,k) &
2389  + b_stag_y(j,k+1 ) )
2390 
2391  IF ( dabs(b_cent(j,k)-b_avg) .GT. 1.d-5 ) THEN
2392 
2393  WRITE(*,*) 'j,k',j,k
2394  WRITE(*,*) 'B_cent(j,k),B_avg1',b_cent(j,k),b_avg
2395  WRITE(*,*)
2396  READ(*,*)
2397 
2398  END IF
2399 
2400  ! Forth check for bi-linearity of recontructed B
2401  b_avg = 0.25d0 * ( b_ver(j,k) + b_ver(j+1,k ) + b_ver(j,k+1) &
2402  + b_ver(j+1,k+1) )
2403 
2404  IF ( dabs(b_cent(j,k)-b_avg) .GT. 1.d-5 ) THEN
2405 
2406  WRITE(*,*) 'j,k',j,k
2407  WRITE(*,*) 'B_cent(j,k),B_avg2',b_cent(j,k),b_avg
2408  WRITE(*,*)
2409  READ(*,*)
2410 
2411  END IF
2412 
2413  ! Update the solution with erosion/deposition terms
2414  q(1:n_eqns,j,k) = q(1:n_eqns,j,k) + dt * eqns_term(1:n_eqns)
2415 
2416  IF ( solid_transport_flag ) q(4,j,k) = min(q(4,j,k),q(1,j,k))
2417 
2418  ! Check for negative thickness
2419  IF ( q(1,j,k) .LT. 0.d0 ) THEN
2420 
2421  WRITE(*,*) 'j,k',j,k
2422  WRITE(*,*) 'dt',dt
2423  WRITE(*,*) 'before erosion'
2424  WRITE(*,*) 'qc',q(1:n_eqns,j,k) - dt * eqns_term(1:n_eqns)
2425  WRITE(*,*) 'B_cent',b_cent(j,k) - dt * topo_term
2426  CALL qc_to_qp(q(1:n_vars,j,k) - dt * eqns_term(1:n_eqns), &
2427  b_cent(j,k) - dt*( avg_deposition_term - avg_erosion_term ), &
2428  qp(1:n_vars,j,k))
2429  WRITE(*,*) 'qp',qp(1:n_vars,j,k)
2430 
2431  WRITE(*,*) 'after update erosion/deposition'
2432  WRITE(*,*) vertex_erosion_terms(j,k) , &
2433  vertex_erosion_terms(j+1,k) , &
2434  vertex_erosion_terms(j,k+1) , &
2435  vertex_erosion_terms(j+1,k+1)
2436 
2437  WRITE(*,*) 'deposition at vertexes', vertex_deposition_terms(j,k) ,&
2438  vertex_deposition_terms(j+1,k) , &
2439  vertex_deposition_terms(j,k+1) , &
2440  vertex_deposition_terms(j+1,k+1)
2441 
2442  WRITE(*,*)
2443 
2444  WRITE(*,*) 'avg_deposition_term',avg_deposition_term
2445  WRITE(*,*) 'avg_erosion_term',avg_erosion_term
2446  WRITE(*,*) 'qc',q(1:n_vars,j,k)
2447 
2448  READ(*,*)
2449 
2450  END IF
2451 
2452 
2453 
2454  IF ( verbose_level .GE. 2 ) THEN
2455 
2456  WRITE(*,*) 'avg_deposition_term , avg_erosion_term', &
2457  avg_deposition_term , avg_erosion_term
2458 
2459  WRITE(*,*) 'after update erosion/deposition: j,k,q(:,j,k),B(j,k)', &
2460  j,k,q(:,j,k),b_cent(j,k)
2461 
2462  READ(*,*)
2463 
2464  END IF
2465 
2466  END DO
2467 
2468  END DO
2469 
2470  RETURN
2471 
2472  END SUBROUTINE update_erosion_deposition_ver
2473 
2474 
2475  !******************************************************************************
2477  !
2480  !
2483  !
2487  !******************************************************************************
2488 
2489  SUBROUTINE eval_explicit_terms( q_expl , expl_terms )
2490 
2491  USE constitutive_2d, ONLY : eval_expl_terms
2492 
2493  IMPLICIT NONE
2494 
2495  REAL*8, INTENT(IN) :: q_expl(n_vars,comp_cells_x,comp_cells_y)
2496  REAL*8, INTENT(OUT) :: expl_terms(n_eqns,comp_cells_x,comp_cells_y)
2497 
2498  REAL*8 :: qc(n_vars)
2499  REAL*8 :: expl_forces_term(n_eqns)
2500 
2501  INTEGER :: j,k
2502 
2503  DO k = 1,comp_cells_y
2504 
2505  DO j = 1,comp_cells_x
2506 
2507  qc = q_expl(1:n_vars,j,k)
2508 
2509  CALL eval_expl_terms( b_prime_x(j,k), b_prime_y(j,k), source_xy(j,k) ,&
2510  qc, expl_forces_term )
2511 
2512  expl_terms(1:n_eqns,j,k) = expl_forces_term
2513 
2514  ENDDO
2515 
2516  END DO
2517 
2518  END SUBROUTINE eval_explicit_terms
2519 
2520  !******************************************************************************
2522  !
2527  !
2530  !
2534  !******************************************************************************
2535 
2536  SUBROUTINE eval_hyperbolic_terms( q_expl , divFlux )
2537 
2538  ! External variables
2539  USE geometry_2d, ONLY : dx,dy
2540  USE parameters_2d, ONLY : solver_scheme
2541 
2542  IMPLICIT NONE
2543 
2544  REAL*8, INTENT(IN) :: q_expl(n_vars,comp_cells_x,comp_cells_y)
2545  REAL*8, INTENT(OUT) :: divFlux(n_eqns,comp_cells_x,comp_cells_y)
2546 
2547  REAL*8 :: q_old(n_vars,comp_cells_x,comp_cells_y)
2548 
2549  REAL*8 :: h_new , h_old
2550 
2551  REAL*8 :: tcpu0,tcpu1,tcpu2,tcpu3,tcpu4
2552 
2553  INTEGER :: i, j, k
2554 
2555  q_old = q
2556 
2557  q = q_expl
2558 
2559  CALL cpu_time(tcpu0)
2560 
2561  ! Linear reconstruction of the physical variables at the interfaces
2562  CALL reconstruction
2563 
2564  CALL cpu_time(tcpu1)
2565  !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu1-tcpu0,'seconds'
2566 
2567 
2568  ! Evaluation of the maximum local speeds at the interfaces
2569  CALL eval_speeds
2570 
2571  CALL cpu_time(tcpu2)
2572  !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu2-tcpu1,'seconds'
2573 
2574  ! Evaluation of the numerical fluxes
2575  SELECT CASE ( solver_scheme )
2576 
2577  CASE ("LxF")
2578 
2579  CALL eval_flux_lxf
2580 
2581  CASE ("GFORCE")
2582 
2583  CALL eval_flux_gforce
2584 
2585  CASE ("KT")
2586 
2587  CALL eval_flux_kt
2588 
2589  END SELECT
2590 
2591  CALL cpu_time(tcpu3)
2592  !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu3-tcpu2,'seconds'
2593 
2594  ! Advance in time the solution
2595  DO k = 1,comp_cells_y
2596 
2597  DO j = 1,comp_cells_x
2598 
2599  DO i=1,n_eqns
2600 
2601  divflux(i,j,k) = 0.d0
2602 
2603  IF ( comp_cells_x .GT. 1 ) THEN
2604 
2605  divflux(i,j,k) = divflux(i,j,k) + &
2606  ( h_interface_x(i,j+1,k) - h_interface_x(i,j,k) ) / dx
2607 
2608  END IF
2609 
2610  IF ( comp_cells_y .GT. 1 ) THEN
2611 
2612  divflux(i,j,k) = divflux(i,j,k) + &
2613  ( h_interface_y(i,j,k+1) - h_interface_y(i,j,k) ) / dy
2614 
2615  END IF
2616 
2617  END DO
2618 
2619  h_old = q_expl(1,j,k)
2620  h_new = h_old - dt * divflux(1,j,k)
2621 
2622  ENDDO
2623 
2624  END DO
2625 
2626  CALL cpu_time(tcpu4)
2627  !WRITE(*,*) 'eval_hyperbolic_terms: Time taken by the code was',tcpu4-tcpu4,'seconds'
2628 
2629  q = q_old
2630 
2631  END SUBROUTINE eval_hyperbolic_terms
2632 
2633 
2634  !******************************************************************************
2636  !
2642  !******************************************************************************
2643 
2644  SUBROUTINE eval_flux_kt
2645 
2646  ! External procedures
2647  USE constitutive_2d, ONLY : eval_fluxes
2648 
2649  IMPLICIT NONE
2650 
2651  REAL*8 :: fluxL(n_eqns)
2652  REAL*8 :: fluxR(n_eqns)
2653  REAL*8 :: fluxB(n_eqns)
2654  REAL*8 :: fluxT(n_eqns)
2655 
2656  REAL*8 :: flux_avg_x(n_eqns)
2657  REAL*8 :: flux_avg_y(n_eqns)
2658 
2659  INTEGER :: i,j,k
2660 
2661  h_interface_x = 0.d0
2662  h_interface_y = 0.d0
2663 
2664  IF ( comp_cells_x .GT. 1 ) THEN
2665 
2666  DO k = 1,comp_cells_y
2667 
2668  DO j = 1,comp_interfaces_x
2669 
2670  CALL eval_fluxes( r_qj = q_interfacel(1:n_vars,j,k) , &
2671  r_flux=fluxl , dir=1 )
2672 
2673  CALL eval_fluxes( r_qj = q_interfacer(1:n_vars,j,k) , &
2674  r_flux=fluxr , dir=1 )
2675 
2676  CALL average_kt( a_interface_xneg(:,j,k), a_interface_xpos(:,j,k) ,&
2677  fluxl , fluxr , flux_avg_x )
2678 
2679  eqns_loop:DO i=1,n_eqns
2680 
2681  IF ( a_interface_xneg(i,j,k) .EQ. a_interface_xpos(i,j,k) ) THEN
2682 
2683  h_interface_x(i,j,k) = 0.d0
2684 
2685  ELSE
2686 
2687  h_interface_x(i,j,k) = flux_avg_x(i) &
2688  + ( a_interface_xpos(i,j,k) * a_interface_xneg(i,j,k) ) &
2689  / ( a_interface_xpos(i,j,k) - a_interface_xneg(i,j,k) ) &
2690  * ( q_interfacer(i,j,k) - q_interfacel(i,j,k) )
2691 
2692  END IF
2693 
2694  ENDDO eqns_loop
2695 
2696  ! In the equation for mass and for trasnport (T,alphas) if the
2697  ! velocities at the interfaces are null, then the flux is null
2698  IF ( ( q_interfacel(2,j,k) .EQ. 0.d0 ) .AND. &
2699  ( q_interfacer(2,j,k) .EQ. 0.d0 ) ) THEN
2700 
2701  h_interface_x(1,j,k) = 0.d0
2702  h_interface_x(4:n_vars,j,k) = 0.d0
2703 
2704  END IF
2705 
2706  END DO
2707 
2708  END DO
2709 
2710  END IF
2711 
2712  IF ( comp_cells_y .GT. 1 ) THEN
2713 
2714  DO k = 1,comp_interfaces_y
2715 
2716  DO j = 1,comp_cells_x
2717 
2718  CALL eval_fluxes( r_qj = q_interfaceb(1:n_vars,j,k) , &
2719  r_flux=fluxb , dir=2 )
2720 
2721  CALL eval_fluxes( r_qj = q_interfacet(1:n_vars,j,k) , &
2722  r_flux=fluxt , dir=2 )
2723 
2724  CALL average_kt( a_interface_yneg(:,j,k) , &
2725  a_interface_ypos(:,j,k) , fluxb , fluxt , flux_avg_y )
2726 
2727  DO i=1,n_eqns
2728 
2729  IF ( a_interface_yneg(i,j,k) .EQ. a_interface_ypos(i,j,k) ) THEN
2730 
2731  h_interface_y(i,j,k) = 0.d0
2732 
2733  ELSE
2734 
2735  h_interface_y(i,j,k) = flux_avg_y(i) &
2736  + ( a_interface_ypos(i,j,k) * a_interface_yneg(i,j,k) ) &
2737  / ( a_interface_ypos(i,j,k) - a_interface_yneg(i,j,k) ) &
2738  * ( q_interfacet(i,j,k) - q_interfaceb(i,j,k) )
2739 
2740  END IF
2741 
2742  END DO
2743 
2744  ! In the equation for mass and for trasnport (T,alphas) if the
2745  ! velocities at the interfaces are null, then the flux is null
2746  IF ( ( q_interfaceb(3,j,k) .EQ. 0.d0 ) .AND. &
2747  ( q_interfacet(3,j,k) .EQ. 0.d0 ) ) THEN
2748 
2749  h_interface_y(1,j,k) = 0.d0
2750  h_interface_y(4:n_vars,j,k) = 0.d0
2751 
2752  END IF
2753 
2754  ENDDO
2755 
2756  END DO
2757 
2758  END IF
2759 
2760  END SUBROUTINE eval_flux_kt
2761 
2762  !******************************************************************************
2764  !
2775  !******************************************************************************
2776 
2777 
2778  SUBROUTINE average_kt( a1 , a2 , w1 , w2 , w_avg )
2779 
2780  IMPLICIT NONE
2781 
2782  REAL*8, INTENT(IN) :: a1(:) , a2(:)
2783  REAL*8, INTENT(IN) :: w1(:) , w2(:)
2784  REAL*8, INTENT(OUT) :: w_avg(:)
2785 
2786  INTEGER :: n
2787  INTEGER :: i
2788 
2789  n = SIZE( a1 )
2790 
2791  DO i=1,n
2792 
2793  IF ( a1(i) .EQ. a2(i) ) THEN
2794 
2795  w_avg(i) = 0.5d0 * ( w1(i) + w2(i) )
2796  w_avg(i) = 0.d0
2797 
2798  ELSE
2799 
2800  w_avg(i) = ( a2(i) * w1(i) - a1(i) * w2(i) ) / ( a2(i) - a1(i) )
2801 
2802  END IF
2803 
2804  END DO
2805 
2806  END SUBROUTINE average_kt
2807 
2808  !******************************************************************************
2813  !******************************************************************************
2814 
2815  SUBROUTINE eval_flux_gforce
2816 
2817  ! to be implemented
2818  WRITE(*,*) 'method not yet implemented in 2-d case'
2819 
2820  END SUBROUTINE eval_flux_gforce
2821 
2822  !******************************************************************************
2827  !******************************************************************************
2828 
2829  SUBROUTINE eval_flux_lxf
2830 
2831  ! to be implemented
2832  WRITE(*,*) 'method not yet implemented in 2-d case'
2833 
2834  END SUBROUTINE eval_flux_lxf
2835 
2836 
2837  !******************************************************************************
2839  !
2846  !******************************************************************************
2847 
2848  SUBROUTINE reconstruction
2849 
2850  ! External procedures
2851  USE constitutive_2d, ONLY : qc_to_qp , qp_to_qc
2852  USE constitutive_2d, ONLY : qc_to_qc2 , qc2_to_qc
2853  USE parameters_2d, ONLY : limiter
2854 
2855  ! External variables
2856  USE geometry_2d, ONLY : x_comp , x_stag , y_comp , y_stag , dx , dx2 , dy , &
2857  dy2
2858 
2859  USE parameters_2d, ONLY : reconstr_variables , reconstr_coeff
2860 
2861  IMPLICIT NONE
2862 
2863  REAL*8 :: qc(n_vars)
2864  REAL*8 :: qrec(n_vars,comp_cells_x,comp_cells_y)
2865  REAL*8 :: qrecW(n_vars)
2866  REAL*8 :: qrecE(n_vars)
2867  REAL*8 :: qrecS(n_vars)
2868  REAL*8 :: qrecN(n_vars)
2869  REAL*8 :: qrec_bdry(n_vars)
2870 
2871  REAL*8 :: qrec_stencil(3)
2872  REAL*8 :: x_stencil(3)
2873  REAL*8 :: y_stencil(3)
2874  REAL*8 :: qrec_prime_x
2875  REAL*8 :: qrec_prime_y
2876 
2877  INTEGER :: j,k
2878  INTEGER :: i
2879 
2880  ! Compute the variable to reconstruct (phys or cons)
2881  DO k = 1,comp_cells_y
2882 
2883  DO j = 1,comp_cells_x
2884 
2885  qc = q(1:n_vars,j,k)
2886 
2887  IF ( reconstr_variables .EQ. 'phys' ) THEN
2888 
2889  CALL qc_to_qp( qc , b_cent(j,k) , qrec(1:n_vars,j,k) )
2890 
2891  ELSEIF ( reconstr_variables .EQ. 'cons' ) THEN
2892 
2893  CALL qc_to_qc2( qc , b_cent(j,k) , qrec(1:n_vars,j,k) )
2894 
2895  END IF
2896 
2897  IF ( solid_transport_flag ) THEN
2898 
2899  IF ( qrec(4,j,k) .GT. 1.d0 ) THEN
2900 
2901  WRITE(*,*) 'reconstruction: j,k',j,k
2902  WRITE(*,*) 'qrec(4,j,k)',qrec(4,j,k)
2903  WRITE(*,*) 'q(1:n_vars,j,k)',q(1:n_vars,j,k)
2904  WRITE(*,*) 'B_cent(j,k)', b_cent(j,k)
2905  WRITE(*,*) 'h',q(1,j,k)-b_cent(j,k)
2906  READ(*,*)
2907 
2908  END IF
2909 
2910  END IF
2911 
2912  END DO
2913 
2914  ENDDO
2915 
2916  ! Linear reconstruction
2917 
2918  y_loop:DO k = 1,comp_cells_y
2919 
2920  x_loop:DO j = 1,comp_cells_x
2921 
2922  vars_loop:DO i=1,n_vars
2923 
2924  ! x direction
2925  check_comp_cells_x:IF ( comp_cells_x .GT. 1 ) THEN
2926 
2927  ! west boundary
2928  check_x_boundary:IF (j.EQ.1) THEN
2929 
2930  IF ( bcw(i)%flag .EQ. 0 ) THEN
2931 
2932  x_stencil(1) = x_stag(1)
2933  x_stencil(2:3) = x_comp(1:2)
2934 
2935  qrec_stencil(1) = bcw(i)%value
2936  qrec_stencil(2:3) = qrec(i,1:2,k)
2937 
2938  CALL limit( qrec_stencil , x_stencil , limiter(i) , &
2939  qrec_prime_x )
2940 
2941  ELSEIF ( bcw(i)%flag .EQ. 1 ) THEN
2942 
2943  qrec_prime_x = bcw(i)%value
2944 
2945  ELSEIF ( bcw(i)%flag .EQ. 2 ) THEN
2946 
2947  qrec_prime_x = ( qrec(i,2,k) - qrec(i,1,k) ) / dx
2948 
2949  END IF
2950 
2951  !east boundary
2952  ELSEIF (j.EQ.comp_cells_x) THEN
2953 
2954  IF ( bce(i)%flag .EQ. 0 ) THEN
2955 
2956  qrec_stencil(3) = bce(i)%value
2957  qrec_stencil(1:2) = qrec(i,comp_cells_x-1:comp_cells_x,k)
2958 
2959  x_stencil(3) = x_stag(comp_interfaces_x)
2960  x_stencil(1:2) = x_comp(comp_cells_x-1:comp_cells_x)
2961 
2962  CALL limit( qrec_stencil , x_stencil , limiter(i) , &
2963  qrec_prime_x )
2964 
2965  ELSEIF ( bce(i)%flag .EQ. 1 ) THEN
2966 
2967  qrec_prime_x = bce(i)%value
2968 
2969  ELSEIF ( bce(i)%flag .EQ. 2 ) THEN
2970 
2971  qrec_prime_x = ( qrec(i,comp_cells_x,k) - &
2972  qrec(i,comp_cells_x-1,k) ) / dx
2973 
2974  END IF
2975 
2976  ! internal x cells
2977  ELSE
2978 
2979  x_stencil(1:3) = x_comp(j-1:j+1)
2980  qrec_stencil = qrec(i,j-1:j+1,k)
2981 
2982  CALL limit( qrec_stencil , x_stencil , limiter(i) , &
2983  qrec_prime_x )
2984 
2985  ENDIF check_x_boundary
2986 
2987  qrecw(i) = qrec(i,j,k) - reconstr_coeff * dx2 * qrec_prime_x
2988  qrece(i) = qrec(i,j,k) + reconstr_coeff * dx2 * qrec_prime_x
2989 
2990  ! positivity preserving reconstruction for h (i=1, 1st variable)
2991  IF ( i .EQ. 1 ) THEN
2992 
2993  IF ( qrece(i) .LT. b_stag_x(j+1,k) ) THEN
2994 
2995  qrec_prime_x = ( b_stag_x(j+1,k) - qrec(i,j,k) ) / dx2
2996 
2997  qrece(i) = qrec(i,j,k) + dx2 * qrec_prime_x
2998 
2999  qrecw(i) = 2.d0 * qrec(i,j,k) - qrece(i)
3000 
3001  ENDIF
3002 
3003  IF ( qrecw(i) .LT. b_stag_x(j,k) ) THEN
3004 
3005  qrec_prime_x = ( qrec(i,j,k) - b_stag_x(j,k) ) / dx2
3006 
3007  qrecw(i) = qrec(i,j,k) - dx2 * qrec_prime_x
3008 
3009  qrece(i) = 2.d0 * qrec(i,j,k) - qrecw(i)
3010 
3011  ENDIF
3012 
3013  END IF
3014 
3015  END IF check_comp_cells_x
3016 
3017  ! y-direction
3018  check_comp_cells_y:IF ( comp_cells_y .GT. 1 ) THEN
3019 
3020  ! South boundary
3021  check_y_boundary:IF (k.EQ.1) THEN
3022 
3023  IF ( bcs(i)%flag .EQ. 0 ) THEN
3024 
3025  qrec_stencil(1) = bcs(i)%value
3026  qrec_stencil(2:3) = qrec(i,j,1:2)
3027 
3028  y_stencil(1) = y_stag(1)
3029  y_stencil(2:3) = y_comp(1:2)
3030 
3031  CALL limit( qrec_stencil , y_stencil , limiter(i) , &
3032  qrec_prime_y )
3033 
3034  ELSEIF ( bcs(i)%flag .EQ. 1 ) THEN
3035 
3036  qrec_prime_y = bcs(i)%value
3037 
3038  ELSEIF ( bcs(i)%flag .EQ. 2 ) THEN
3039 
3040  qrec_prime_y = ( qrec(i,j,2) - qrec(i,j,1) ) / dy
3041 
3042  END IF
3043 
3044  ! North boundary
3045  ELSEIF ( k .EQ. comp_cells_y ) THEN
3046 
3047  IF ( bcn(i)%flag .EQ. 0 ) THEN
3048 
3049  qrec_stencil(3) = bcn(i)%value
3050  qrec_stencil(1:2) = qrec(i,j,comp_cells_y-1:comp_cells_y)
3051 
3052  y_stencil(3) = y_stag(comp_interfaces_y)
3053  y_stencil(1:2) = y_comp(comp_cells_y-1:comp_cells_y)
3054 
3055  CALL limit( qrec_stencil , y_stencil , limiter(i) , &
3056  qrec_prime_y )
3057 
3058  ELSEIF ( bcn(i)%flag .EQ. 1 ) THEN
3059 
3060  qrec_prime_y = bcn(i)%value
3061 
3062  ELSEIF ( bcn(i)%flag .EQ. 2 ) THEN
3063 
3064  qrec_prime_y = ( qrec(i,j,comp_cells_y) - &
3065  qrec(i,j,comp_cells_y-1) ) / dy
3066 
3067  END IF
3068 
3069  ! Internal y cells
3070  ELSE
3071 
3072  y_stencil(1:3) = y_comp(k-1:k+1)
3073  qrec_stencil = qrec(i,j,k-1:k+1)
3074 
3075  CALL limit( qrec_stencil , y_stencil , limiter(i) , &
3076  qrec_prime_y )
3077 
3078  ENDIF check_y_boundary
3079 
3080  qrecs(i) = qrec(i,j,k) - reconstr_coeff * dy2 * qrec_prime_y
3081  qrecn(i) = qrec(i,j,k) + reconstr_coeff * dy2 * qrec_prime_y
3082 
3083  ! positivity preserving reconstruction for h
3084  IF ( i .EQ. 1 ) THEN
3085 
3086  IF ( qrecn(i) .LT. b_stag_y(j,k+1) ) THEN
3087 
3088  qrec_prime_y = ( b_stag_y(j,k+1) - qrec(i,j,k) ) / dy2
3089 
3090  qrecn(i) = qrec(i,j,k) + dy2 * qrec_prime_y
3091 
3092  qrecs(i) = 2.d0 * qrec(i,j,k) - qrecn(i)
3093 
3094  ENDIF
3095 
3096  IF ( qrecs(i) .LT. b_stag_y(j,k) ) THEN
3097 
3098  qrec_prime_y = ( qrec(i,j,k) - b_stag_y(j,k) ) / dy2
3099 
3100  qrecs(i) = qrec(i,j,k) - dy2 * qrec_prime_y
3101 
3102  qrecn(i) = 2.d0 * qrec(i,j,k) - qrecs(i)
3103 
3104  ENDIF
3105 
3106  END IF
3107 
3108  ENDIF check_comp_cells_y
3109 
3110  ENDDO vars_loop
3111 
3112 
3113  IF ( comp_cells_x .GT. 1 ) THEN
3114 
3115  IF ( reconstr_variables .EQ. 'phys' ) THEN
3116 
3117  CALL qp_to_qc( qrecw , b_stag_x(j,k) , q_interfacer(:,j,k) )
3118  CALL qp_to_qc( qrece , b_stag_x(j+1,k) , q_interfacel(:,j+1,k) )
3119 
3120  ELSE IF ( reconstr_variables .EQ. 'cons' ) THEN
3121 
3122  CALL qc2_to_qc( qrecw , b_stag_x(j,k), q_interfacer(:,j,k) )
3123  CALL qc2_to_qc( qrece , b_stag_x(j+1,k) , q_interfacel(:,j+1,k) )
3124 
3125  END IF
3126 
3127  IF ( ( j.EQ. -1000 ) .OR. ( j.EQ. -1001 ) )THEN
3128 
3129  WRITE(*,*) 'j',j
3130  WRITE(*,*) 'B_stag_x(j,k)', b_stag_x(j,k)
3131  WRITE(*,*) 'B_cent(j,k)',b_cent(j,k)
3132  WRITE(*,*) 'B_stag_x(j+1,k)', b_stag_x(j+1,k)
3133  WRITE(*,*)
3134  WRITE(*,*) 'qrecW', qrecw
3135  WRITE(*,*) 'qrec', qrec(:,j,k)
3136  WRITE(*,*) 'qrecE', qrece
3137  WRITE(*,*)
3138  WRITE(*,*) 'q_interfaceR(:,j,k)', q_interfacer(:,j,k)
3139  WRITE(*,*) 'q(:,j,k)',q(:,j,k)
3140  WRITE(*,*) 'q_interfaceL(:,j+1,k)',q_interfacel(:,j+1,k)
3141 
3142  READ(*,*)
3143 
3144  END IF
3145 
3146  IF ( q(1,j,k) .EQ. 0.d0 ) THEN
3147 
3148  q_interfacer(1,j,k) = 0.d0
3149  q_interfacel(1,j+1,k) = 0.d0
3150 
3151  IF ( solid_transport_flag ) THEN
3152 
3153  q_interfacer(4,j,k) = 0.d0
3154  q_interfacel(4,j+1,k) = 0.d0
3155 
3156  END IF
3157 
3158  END IF
3159 
3160  ELSE
3161 
3162  q_interfacer(:,j,k) = q(:,j,k)
3163  q_interfacel(:,j+1,k) = q(:,j,k)
3164 
3165  END IF
3166 
3167 
3168  IF ( comp_cells_y .GT. 1 ) THEN
3169 
3170  IF ( reconstr_variables .EQ. 'phys' ) THEN
3171 
3172  CALL qp_to_qc( qrecs , b_stag_y(j,k) , q_interfacet(:,j,k) )
3173  CALL qp_to_qc( qrecn , b_stag_y(j,k+1) , q_interfaceb(:,j,k+1) )
3174 
3175  ELSE IF ( reconstr_variables .EQ. 'cons' ) THEN
3176 
3177  CALL qc2_to_qc( qrecs , b_stag_y(j,k), q_interfacet(:,j,k) )
3178  CALL qc2_to_qc( qrecn , b_stag_y(j,k+1), q_interfaceb(:,j,k+1) )
3179 
3180  END IF
3181 
3182  IF ( q(1,j,k) .EQ. 0.d0 ) THEN
3183 
3184  q_interfacet(1,j,k) = 0.d0
3185  q_interfaceb(1,j,k+1) = 0.d0
3186 
3187  IF ( solid_transport_flag) THEN
3188 
3189  q_interfacet(4,j,k) = 0.d0
3190  q_interfaceb(4,j,k+1) = 0.d0
3191 
3192  END IF
3193 
3194  END IF
3195 
3196  ELSE
3197 
3198  q_interfacet(:,j,k) = q(:,j,k)
3199  q_interfaceb(:,j,k+1) = q(:,j,k)
3200 
3201  END IF
3202 
3203 
3204  ! Boundary conditions outside the domain (external interfaces)
3205  check_if_solve_along_y:IF ( comp_cells_y .GT. 1 ) THEN
3206 
3207  ! South q_interfaceB(:,j,1)
3208  south_bc:IF ( k .EQ. 1 ) THEN
3209 
3210  DO i=1,n_vars
3211 
3212  IF ( bcs(i)%flag .EQ. 0 ) THEN
3213 
3214  qrec_bdry(i) = bcs(i)%value
3215 
3216  ELSE
3217 
3218  qrec_bdry(i) = qrecs(i)
3219 
3220  END IF
3221 
3222  ENDDO
3223 
3224  IF ( reconstr_variables .EQ. 'phys' ) THEN
3225 
3226  CALL qp_to_qc( qrec_bdry, b_stag_y(j,1), q_interfaceb(:,j,1) )
3227 
3228  ELSEIF ( reconstr_variables .EQ. 'cons' ) THEN
3229 
3230  CALL qc2_to_qc(qrec_bdry, b_stag_y(j,1), q_interfaceb(:,j,1) )
3231 
3232  END IF
3233 
3234  ENDIF south_bc
3235 
3236  ! North qT(i,j,comp_interfaces_y)
3237  north_bc:IF ( k .EQ. comp_cells_y ) THEN
3238 
3239  DO i=1,n_vars
3240 
3241  IF ( bcn(i)%flag .EQ. 0 ) THEN
3242 
3243  qrec_bdry(i) = bcn(i)%value
3244 
3245  ELSE
3246 
3247  qrec_bdry(i) = qrecn(i)
3248 
3249  END IF
3250 
3251  ENDDO
3252 
3253  IF ( reconstr_variables .EQ. 'phys' ) THEN
3254 
3255  CALL qp_to_qc( qrec_bdry , b_stag_y(j,comp_interfaces_y) , &
3256  q_interfacet(:,j,comp_interfaces_y) )
3257 
3258  ELSEIF ( reconstr_variables .EQ. 'cons' ) THEN
3259 
3260  CALL qc2_to_qc( qrec_bdry , b_stag_y(j,comp_interfaces_y) , &
3261  q_interfacet(:,j,comp_interfaces_y) )
3262 
3263  END IF
3264 
3265  ENDIF north_bc
3266 
3267  ELSE
3268  ! 1D-simulation in the x-direction
3269 
3270  IF ( k .EQ. 1 ) THEN
3271 
3272  q_interfaceb(:,j,1) = q(:,j,k)
3273 
3274  END IF
3275 
3276  IF ( k .EQ. comp_cells_y ) THEN
3277 
3278  q_interfacet(:,j,comp_interfaces_y) = q(:,j,k)
3279 
3280  END IF
3281 
3282  END IF check_if_solve_along_y
3283 
3284  check_if_solve_along_x:IF ( comp_cells_x .GT. 1 ) THEN
3285 
3286  ! West q_interfaceL(:,1,k)
3287  west_bc:IF ( j.EQ.1 ) THEN
3288 
3289  DO i=1,n_vars
3290 
3291  IF ( bcw(i)%flag .EQ. 0 ) THEN
3292 
3293  qrec_bdry(i) = bcw(i)%value
3294 
3295  ELSE
3296 
3297  qrec_bdry(i) = qrecw(i)
3298 
3299  END IF
3300 
3301  ENDDO
3302 
3303  IF ( reconstr_variables .EQ. 'phys' ) THEN
3304 
3305  CALL qp_to_qc( qrec_bdry, b_stag_x(1,k), q_interfacel(:,1,k) )
3306 
3307  ELSEIF ( reconstr_variables .EQ. 'cons' ) THEN
3308 
3309  CALL qc2_to_qc( qrec_bdry, b_stag_x(1,k), q_interfacel(:,1,k) )
3310 
3311  END IF
3312 
3313  q_interfacer(:,1,k) = q_interfacel(:,1,k)
3314 
3315  ENDIF west_bc
3316 
3317  ! East q_interfaceR(:,comp_interfaces_x,k)
3318  east_bc:IF ( j.EQ.comp_cells_x ) THEN
3319 
3320  DO i=1,n_vars
3321 
3322  IF ( bce(i)%flag .EQ. 0 ) THEN
3323 
3324  qrec_bdry(i) = bce(i)%value
3325 
3326  ELSE
3327 
3328  qrec_bdry(i) = qrece(i)
3329 
3330  END IF
3331 
3332  ENDDO
3333 
3334  IF ( reconstr_variables .EQ. 'phys' ) THEN
3335 
3336  CALL qp_to_qc( qrec_bdry , b_stag_x(comp_interfaces_x,k) , &
3337  q_interfacer(:,comp_interfaces_x,k) )
3338 
3339  ELSEIF ( reconstr_variables .EQ. 'cons' ) THEN
3340 
3341  CALL qc2_to_qc( qrec_bdry , b_stag_x(comp_interfaces_x,k) , &
3342  q_interfacer(:,comp_interfaces_x,k) )
3343 
3344  END IF
3345 
3346  q_interfacel(:,comp_interfaces_x,k) = &
3347  q_interfacer(:,comp_interfaces_x,k)
3348 
3349  ENDIF east_bc
3350 
3351  ELSE
3352 
3353  ! 1D-simulation in the y-direction
3354  IF ( j .EQ. 1 ) THEN
3355 
3356  q_interfacel(:,1,k) = q(:,j,k)
3357 
3358  END IF
3359 
3360  IF ( j .EQ. comp_cells_x ) THEN
3361 
3362  q_interfacer(:,comp_interfaces_x,k) = q(:,j,k)
3363 
3364  END IF
3365 
3366  END IF check_if_solve_along_x
3367 
3368  END DO x_loop
3369 
3370  END DO y_loop
3371 
3372  END SUBROUTINE reconstruction
3373 
3374 
3375  !******************************************************************************
3377  !
3383  !******************************************************************************
3384 
3385  SUBROUTINE eval_speeds
3386 
3387  ! External procedures
3388  USE constitutive_2d, ONLY : eval_local_speeds_x, eval_local_speeds_y
3389 
3390  IMPLICIT NONE
3391 
3392  REAL*8 :: abslambdaL_min(n_vars) , abslambdaL_max(n_vars)
3393  REAL*8 :: abslambdaR_min(n_vars) , abslambdaR_max(n_vars)
3394  REAL*8 :: abslambdaB_min(n_vars) , abslambdaB_max(n_vars)
3395  REAL*8 :: abslambdaT_min(n_vars) , abslambdaT_max(n_vars)
3396  REAL*8 :: min_r(n_vars) , max_r(n_vars)
3397 
3398  INTEGER :: j,k
3399 
3400  IF ( comp_cells_x .GT. 1 ) THEN
3401 
3402  DO j = 1,comp_interfaces_x
3403 
3404  DO k = 1, comp_cells_y
3405 
3406  CALL eval_local_speeds_x( q_interfacel(:,j,k) , abslambdal_min , &
3407  abslambdal_max )
3408 
3409  CALL eval_local_speeds_x( q_interfacer(:,j,k) , abslambdar_min , &
3410  abslambdar_max )
3411 
3412  min_r = min(abslambdal_min , abslambdar_min , 0.0d0)
3413  max_r = max(abslambdal_max , abslambdar_max , 0.0d0)
3414 
3415  a_interface_xneg(:,j,k) = min_r
3416  a_interface_xpos(:,j,k) = max_r
3417 
3418  ENDDO
3419 
3420  END DO
3421 
3422  END IF
3423 
3424  IF ( comp_cells_y .GT. 1 ) THEN
3425 
3426  DO j = 1,comp_cells_x
3427 
3428  DO k = 1,comp_interfaces_y
3429 
3430  CALL eval_local_speeds_y( q_interfaceb(:,j,k) , abslambdab_min , &
3431  abslambdab_max )
3432 
3433  CALL eval_local_speeds_y( q_interfacet(:,j,k) , abslambdat_min , &
3434  abslambdat_max )
3435 
3436  min_r = min(abslambdab_min , abslambdat_min , 0.0d0)
3437  max_r = max(abslambdab_max , abslambdat_max , 0.0d0)
3438 
3439  a_interface_yneg(:,j,k) = min_r
3440  a_interface_ypos(:,j,k) = max_r
3441 
3442  ENDDO
3443 
3444  END DO
3445 
3446  END IF
3447 
3448  END SUBROUTINE eval_speeds
3449 
3450 
3451  !******************************************************************************
3453  !
3468  !******************************************************************************
3469 
3470  SUBROUTINE limit( v , z , limiter , slope_lim )
3471 
3472  USE parameters_2d, ONLY : theta
3473 
3474  IMPLICIT none
3475 
3476  REAL*8, INTENT(IN) :: v(3)
3477  REAL*8, INTENT(IN) :: z(3)
3478  INTEGER, INTENT(IN) :: limiter
3479 
3480  REAL*8, INTENT(OUT) :: slope_lim
3481 
3482  REAL*8 :: a , b , c
3483 
3484  REAL*8 :: sigma1 , sigma2
3485 
3486  a = ( v(3) - v(2) ) / ( z(3) - z(2) )
3487  b = ( v(2) - v(1) ) / ( z(2) - z(1) )
3488  c = ( v(3) - v(1) ) / ( z(3) - z(1) )
3489 
3490  SELECT CASE (limiter)
3491 
3492  CASE ( 0 )
3493 
3494  slope_lim = 0.d0
3495 
3496  CASE ( 1 )
3497 
3498  ! minmod
3499 
3500  slope_lim = minmod(a,b)
3501 
3502  CASE ( 2 )
3503 
3504  ! superbee
3505 
3506  sigma1 = minmod( a , 2.d0*b )
3507  sigma2 = minmod( 2.d0*a , b )
3508  slope_lim = maxmod( sigma1 , sigma2 )
3509 
3510  CASE ( 3 )
3511 
3512  ! generalized minmod
3513 
3514  slope_lim = minmod( c , theta * minmod( a , b ) )
3515 
3516  CASE ( 4 )
3517 
3518  ! monotonized central-difference (MC, LeVeque p.112)
3519 
3520  slope_lim = minmod( c , 2.0 * minmod( a , b ) )
3521 
3522  END SELECT
3523 
3524  END SUBROUTINE limit
3525 
3526 
3527  REAL*8 FUNCTION minmod(a,b)
3528 
3529  IMPLICIT none
3530 
3531  REAL*8 :: a , b , sa , sb
3532 
3533  IF ( a*b .EQ. 0.d0 ) THEN
3534 
3535  minmod = 0.d0
3536 
3537  ELSE
3538 
3539  sa = a / abs(a)
3540  sb = b / abs(b)
3541 
3542  minmod = 0.5d0 * ( sa+sb ) * min( abs(a) , abs(b) )
3543 
3544  END IF
3545 
3546  END FUNCTION minmod
3547 
3548  REAL*8 function maxmod(a,b)
3549 
3550  IMPLICIT none
3551 
3552  REAL*8 :: a , b , sa , sb
3553 
3554  IF ( a*b .EQ. 0.d0 ) THEN
3555 
3556  maxmod = 0.d0
3557 
3558  ELSE
3559 
3560  sa = a / abs(a)
3561  sb = b / abs(b)
3562 
3563  maxmod = 0.5d0 * ( sa+sb ) * max( abs(a) , abs(b) )
3564 
3565  END IF
3566 
3567  END function maxmod
3568 
3569 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:50
solver_2d::eval_explicit_terms
subroutine eval_explicit_terms(q_expl, expl_terms)
Evaluate the explicit terms.
Definition: solver_2d.f90:2490
solver_2d::q_cellne
real *8, dimension(:,:,:), allocatable q_cellne
Reconstructed value at the NE corner of cell.
Definition: solver_2d.f90:53
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:3386
constitutive_2d::qc_to_qp
subroutine qc_to_qp(qc, B, qp)
Conservative to physical variables.
Definition: constitutive_2d.f90:323
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:104
geometry_2d::dy
real *8 dy
Control volumes size.
Definition: geometry_2d.f90:71
solver_2d::nhj
real *8, dimension(:,:), allocatable nhj
Local Intermediate non-hyperbolic terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:128
parameters_2d::n_rk
integer n_rk
Runge-Kutta order.
Definition: parameters_2d.f90:125
parameters_2d::tol_rel
real *8, parameter tol_rel
Definition: parameters_2d.f90:130
solver_2d::q_rk
real *8, dimension(:,:,:,:), allocatable q_rk
Intermediate solutions of the Runge-Kutta scheme.
Definition: solver_2d.f90:110
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:37
geometry_2d::b_cent
real *8, dimension(:,:), allocatable b_cent
Topography at the centers of the control volumes.
Definition: geometry_2d.f90:47
geometry_2d::comp_cells_x
integer comp_cells_x
Number of control volumes x in the comp. domain.
Definition: geometry_2d.f90:76
solver_2d::eval_flux_gforce
subroutine eval_flux_gforce
Numerical fluxes GFORCE.
Definition: solver_2d.f90:2816
geometry_2d::grav_surf
real *8, dimension(:,:), allocatable grav_surf
gravity vector wrt surface coordinates for each cell
Definition: geometry_2d.f90:59
solver_2d::divflux
real *8, dimension(:,:,:,:), allocatable divflux
Intermediate hyperbolic terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:113
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:85
solver_2d::source_xy
real *8, dimension(:,:), allocatable source_xy
Array defining fraction of cells affected by source term.
Definition: solver_2d.f90:78
constitutive_2d::eval_expl_terms
subroutine eval_expl_terms(Bprimej_x, Bprimej_y, source_xy, qj, expl_term)
Explicit Forces term.
Definition: constitutive_2d.f90:1200
solver_2d::average_kt
subroutine average_kt(a1, a2, w1, w2, w_avg)
averaged KT flux
Definition: solver_2d.f90:2779
parameters_2d
Parameters.
Definition: parameters_2d.f90:7
geometry_2d::dx
real *8 dx
Control volumes size.
Definition: geometry_2d.f90:68
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:119
solver_2d::q_cellse
real *8, dimension(:,:,:), allocatable q_cellse
Reconstructed value at the SE corner of cell.
Definition: solver_2d.f90:57
solver_2d::maxmod
real *8 function maxmod(a, b)
Definition: solver_2d.f90:3549
constitutive_2d::eval_erosion_dep_term
subroutine eval_erosion_dep_term(qj, erosion_term, deposition_term)
Deposition term.
Definition: constitutive_2d.f90:1283
geometry_2d::b_sw
real *8, dimension(:,:), allocatable b_sw
Topography interpolated at the SW corner of the control volumes.
Definition: geometry_2d.f90:38
geometry_2d::comp_cells_y
integer comp_cells_y
Number of control volumes y in the comp. domain.
Definition: geometry_2d.f90:78
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:63
solver_2d::deallocate_solver_variables
subroutine deallocate_solver_variables
Memory deallocation.
Definition: solver_2d.f90:370
solver_2d::expl_terms
real *8, dimension(:,:,:,:), allocatable expl_terms
Intermediate explicit terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:122
parameters_2d::bcw
type(bc), dimension(:), allocatable bcw
bcW&flag defines the west boundary condition:
Definition: parameters_2d.f90:172
solver_2d::nh
real *8, dimension(:,:,:,:), allocatable nh
Intermediate non-hyperbolic terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:116
constitutive_2d::eval_nh_semi_impl_terms
subroutine eval_nh_semi_impl_terms(grav3_surf, c_qj, c_nh_semi_impl_term, r_qj, r_nh_semi_impl_term)
Non-Hyperbolic semi-implicit terms.
Definition: constitutive_2d.f90:1048
solver_2d::solve_mask0
logical, dimension(:,:), allocatable solve_mask0
Definition: solver_2d.f90:80
solver_2d::q_cellsw
real *8, dimension(:,:,:), allocatable q_cellsw
Reconstructed value at the SW corner of cell.
Definition: solver_2d.f90:55
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:101
solver_2d::dt
real *8 dt
Time step.
Definition: solver_2d.f90:83
parameters_2d::n_vars
integer n_vars
Number of conservative variables.
Definition: parameters_2d.f90:120
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:67
solver_2d::eval_flux_lxf
subroutine eval_flux_lxf
Numerical fluxes Lax-Friedrichs.
Definition: solver_2d.f90:2830
geometry_2d::comp_interfaces_y
integer comp_interfaces_y
Number of interfaces (comp_cells_y+1)
Definition: geometry_2d.f90:79
parameters_2d::cfl
real *8 cfl
Courant-Friedrichs-Lewy parameter.
Definition: parameters_2d.f90:18
solver_2d::eval_hyperbolic_terms
subroutine eval_hyperbolic_terms(q_expl, divFlux)
Semidiscrete finite volume central scheme.
Definition: solver_2d.f90:2537
parameters_2d::max_nl_iter
integer, parameter max_nl_iter
Definition: parameters_2d.f90:127
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:146
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:137
solver_2d::update_erosion_deposition_ver
subroutine update_erosion_deposition_ver(dt)
Evaluate the averaged explicit terms.
Definition: solver_2d.f90:1745
solver_2d::mask21
logical, dimension(:,:), allocatable mask21
Definition: solver_2d.f90:85
solver_2d::h_interface_x
real *8, dimension(:,:,:), allocatable h_interface_x
Semidiscrete numerical interface fluxes.
Definition: solver_2d.f90:71
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:92
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:69
constitutive_2d::implicit_flag
logical, dimension(:), allocatable implicit_flag
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:73
solver_2d::minmod
real *8 function minmod(a, b)
Definition: solver_2d.f90:3528
solver_2d::q_interfacel
real *8, dimension(:,:,:), allocatable q_interfacel
Reconstructed value at the left of the x-interface.
Definition: solver_2d.f90:42
solver_2d::omega_tilde
real *8, dimension(:), allocatable omega_tilde
Coefficients for the explicit part of the Runge-Kutta scheme.
Definition: solver_2d.f90:95
parameters_2d::solver_scheme
character(len=20) solver_scheme
Finite volume method: .
Definition: parameters_2d.f90:145
geometry_2d::b_ver
real *8, dimension(:,:), allocatable b_ver
Topography at the vertices of the control volumes.
Definition: geometry_2d.f90:44
constitutive_2d::erosion_coeff
real *8 erosion_coeff
erosion model coefficient
Definition: constitutive_2d.f90:114
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:98
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:131
geometry_2d::b_stag_x
real *8, dimension(:,:), allocatable b_stag_x
Topography at the boundaries (x) of the control volumes.
Definition: geometry_2d.f90:26
solver_2d::q_cellnw
real *8, dimension(:,:,:), allocatable q_cellnw
Reconstructed value at the NW corner of cell.
Definition: solver_2d.f90:51
parameters_2d::theta
real *8 theta
Van Leer limiter parameter.
Definition: parameters_2d.f90:147
solver_2d::timestep2
subroutine timestep2
Time-step computation.
Definition: solver_2d.f90:564
solver_2d::check_solve
subroutine check_solve
Masking of cells to solve.
Definition: solver_2d.f90:443
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:65
constitutive_2d::settling_vel
real *8 settling_vel
hindered settling
Definition: constitutive_2d.f90:111
solver_2d::mask22
logical, dimension(:,:), allocatable mask22
Definition: solver_2d.f90:85
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:1397
geometry_2d::b_ne
real *8, dimension(:,:), allocatable b_ne
Topography interpolated at the NE corner of the control volumes.
Definition: geometry_2d.f90:35
parameters_2d::verbose_level
integer verbose_level
Definition: parameters_2d.f90:155
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:53
solver_2d::timestep
subroutine timestep
Time-step computation.
Definition: solver_2d.f90:487
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:134
geometry_2d::b_stag_y
real *8, dimension(:,:), allocatable b_stag_y
Topography at the boundaries (y) of the control volumes.
Definition: geometry_2d.f90:29
solver_2d::q_interfaceb
real *8, dimension(:,:,:), allocatable q_interfaceb
Reconstructed value at the bottom of the y-interface.
Definition: solver_2d.f90:46
solver_2d::solve_rk_step
subroutine solve_rk_step(Bj, qj, qj_old, a_tilde, a_dirk, a_diag)
Runge-Kutta single step integration.
Definition: solver_2d.f90:1081
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
constitutive_2d::qc2_to_qc
subroutine qc2_to_qc(qc2, B, qc)
Reconstructed to conservative variables.
Definition: constitutive_2d.f90:480
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:90
constitutive_2d::qc_to_qc2
subroutine qc_to_qc2(qc, B, qc2)
Conservative to alternative conservative variables.
Definition: constitutive_2d.f90:451
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:107
constitutive_2d::eval_local_speeds_y
subroutine eval_local_speeds_y(qj, vel_min, vel_max)
Local Characteristic speeds y direction.
Definition: constitutive_2d.f90:291
parameters_2d::reconstr_variables
character(len=20) reconstr_variables
Definition: parameters_2d.f90:22
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:143
solver_2d::q_fv
real *8, dimension(:,:,:), allocatable q_fv
Solution of the finite-volume semidiscrete cheme.
Definition: solver_2d.f90:39
solver_2d::normalize_f
logical normalize_f
Flag for the normalization of the array f in the implicit solution scheme.
Definition: solver_2d.f90:140
parameters_2d::limiter
integer, dimension(10) limiter
Limiter for the slope in the linear reconstruction: .
Definition: parameters_2d.f90:138
parameters_2d::bcn
type(bc), dimension(:), allocatable bcn
bcN&flag defines the north boundary condition:
Definition: parameters_2d.f90:202
geometry_2d::b_se
real *8, dimension(:,:), allocatable b_se
Topography interpolated at the SE corner of the control volumes.
Definition: geometry_2d.f90:41
geometry_2d::dx2
real *8 dx2
Half x Control volumes size.
Definition: geometry_2d.f90:74
geometry_2d::b_nw
real *8, dimension(:,:), allocatable b_nw
Topography interpolated at the NW corner of the control volumes.
Definition: geometry_2d.f90:32
parameters_2d::n_eqns
integer n_eqns
Number of equations.
Definition: parameters_2d.f90:121
solver_2d::eval_jacobian
subroutine eval_jacobian(qj_rel, qj_org, coeff_f, left_matrix)
Evaluate the jacobian.
Definition: solver_2d.f90:1674
constitutive_2d::qp_to_qc
subroutine qp_to_qc(qp, B, qc)
Physical to conservative variables.
Definition: constitutive_2d.f90:372
parameters_2d::tol_abs
real *8, parameter tol_abs
Definition: parameters_2d.f90:129
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:725
constitutive_2d::eval_fluxes
subroutine eval_fluxes(c_qj, r_qj, c_flux, r_flux, dir)
Hyperbolic Fluxes.
Definition: constitutive_2d.f90:563
solver_2d::reconstruction
subroutine reconstruction
Linear reconstruction.
Definition: solver_2d.f90:2849
constitutive_2d::eval_topo_term
subroutine eval_topo_term(qj, deposition_avg_term, erosion_avg_term, eqns_term, topo_term)
Topography modification related term.
Definition: constitutive_2d.f90:1371
solver_2d::t_imex1
real *8 t_imex1
Definition: solver_2d.f90:149
parameters_2d::reconstr_coeff
real *8 reconstr_coeff
Slope coefficient in the linear reconstruction.
Definition: parameters_2d.f90:24
solver_2d::imex_rk_solver
subroutine imex_rk_solver
Runge-Kutta integration.
Definition: solver_2d.f90:634
parameters_2d::bcs
type(bc), dimension(:), allocatable bcs
bcS&flag defines the south boundary condition:
Definition: parameters_2d.f90:192
solver_2d::solve_mask
logical, dimension(:,:), allocatable solve_mask
Definition: solver_2d.f90:80
geometry_2d::comp_interfaces_x
integer comp_interfaces_x
Number of interfaces (comp_cells_x+1)
Definition: geometry_2d.f90:77
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:182
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:1613
solver_2d::q1max
real *8, dimension(:,:), allocatable q1max
Maximum over time of thickness.
Definition: solver_2d.f90:60
solver_2d::limit
subroutine limit(v, z, limiter, slope_lim)
Slope limiter.
Definition: solver_2d.f90:3471
solver_2d::q_interfacer
real *8, dimension(:,:,:), allocatable q_interfacer
Reconstructed value at the right of the x-interface.
Definition: solver_2d.f90:44
solver_2d::divfluxj
real *8, dimension(:,:), allocatable divfluxj
Local Intermediate hyperbolic terms of the Runge-Kutta scheme.
Definition: solver_2d.f90:125
parameters_2d::solid_transport_flag
logical solid_transport_flag
Flag to choose if we model solid phase transport.
Definition: parameters_2d.f90:72
solver_2d::eval_flux_kt
subroutine eval_flux_kt
Semidiscrete numerical fluxes.
Definition: solver_2d.f90:2645
solver_2d::i_rk
integer i_rk
loop counter for the RK iteration
Definition: solver_2d.f90:87
solver_2d::allocate_solver_variables
subroutine allocate_solver_variables
Memory allocation.
Definition: solver_2d.f90:166
solver_2d::t_imex2
real *8 t_imex2
Definition: solver_2d.f90:149
geometry_2d::dy2
real *8 dy2
Half y Control volumes size.
Definition: geometry_2d.f90:75
solver_2d::mask12
logical, dimension(:,:), allocatable mask12
Definition: solver_2d.f90:85
constitutive_2d::eval_local_speeds_x
subroutine eval_local_speeds_x(qj, vel_min, vel_max)
Local Characteristic speeds x direction.
Definition: constitutive_2d.f90:262
solver_2d::q_interfacet
real *8, dimension(:,:,:), allocatable q_interfacet
Reconstructed value at the top of the y-interface.
Definition: solver_2d.f90:48
solver_2d::qp
real *8, dimension(:,:,:), allocatable qp
Physical variables ( )
Definition: solver_2d.f90:75
parameters_2d::n_nh
integer n_nh
Number of non-hyperbolic terms.
Definition: parameters_2d.f90:123
solver_2d::q
real *8, dimension(:,:,:), allocatable q
Conservative variables.
Definition: solver_2d.f90:35
Generated on Thu Apr 4 2019 14:29:02 for LAHARS-MODEL by   doxygen 1.8.11