Changeset 923

Timestamp:

06/19/19 16:06:19 (5 years ago)

Author:

dubos

Message:

devel : really fix bug in tridiagonal solver for vertical diffusion in Venus test case

File:

• codes/icosagcm/devel/src/initial/etat0_venus.f90 (modified) (2 diffs)

codes/icosagcm/devel/src/initial/etat0_venus.f90

@@ -182 +182 @@
     END DO
     ! Back substitution :
-    ! x(i) = D(i)-C(i+1)x(i+1), x(llm)=D(llm)
+    ! x(i) = D(i)-C(i)x(i+1), x(llm)=D(llm)
     !DIR$ SIMD
     DO ij=1,ngrid
@@ -193 +193 @@
        !DIR$ SIMD
        DO ij=1,ngrid
-          xu(ij,l) = xu(ij,l) - C(ij,l+1)*xu(ij,l+1)
-          xv(ij,l) = xv(ij,l) - C(ij,l+1)*xv(ij,l+1)
+          xu(ij,l) = xu(ij,l) - C(ij,l)*xu(ij,l+1)
+          xv(ij,l) = xv(ij,l) - C(ij,l)*xv(ij,l+1)
           du(ij,l) = (xu(ij,l)-u(ij,l))/dt_phys
           dv(ij,l) = (xv(ij,l)-v(ij,l))/dt_phys
        END DO
     END DO
+
+    IF(.FALSE.) THEN
+       ! check correctness of solution of tridiagonal problem
+       PRINT *, 'Max 1 residual of tridiagonal solver', MAXVAL(ABS(Ru))
+       DO l = 2,llm-1
+          !DIR$ SIMD
+          DO ij=1,ngrid
+             Ru(ij,l) = Ru(ij,l) - B(ij,l)*xu(ij,l) + A(ij,l)*xu(ij,l-1) + A(ij,l+1)*xu(ij,l+1)
+          END DO
+       END DO
+       DO ij=1,ngrid
+          Ru(ij,1) = Ru(ij,1) - B(ij,1)*xu(ij,1) + A(ij,2)*xu(ij,2)
+          Ru(ij,llm) = Ru(ij,llm) - B(ij,llm)*xu(ij,llm) + A(ij,llm)*xu(ij,llm-1)
+       END DO
+       PRINT *, 'Max 2 residual of tridiagonal solver', MAXVAL(ABS(Ru))
+    END IF
+
     ! bottom friction + thermal relaxation
     du(:,1)=du(:,1)-kfrict*u(:,1)