Changeset 686 for codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90

Timestamp:

03/08/18 12:48:37 (6 years ago)

Author:

dubos

Message:

devel/unstructured : piecewise-constant vertical remapping

File:

• codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90 (modified) (3 diffs)

codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90

@@ -52 +52 @@
       DO ij=1,primal_num
          ! Compute residual R_il and B_il
-         l=1
          ! bottom interface l=1
-         ml_g2 = gm2*m_il(l,ij)
-         B_il(l,ij) = A_ik(l,ij) + ml_g2 + tau2_g*rho_bot
-         R_il(l,ij) = ml_g2*( Phi_il(l,ij)-Phi_star_il(l,ij)) &
-         + tau2_g*( p_ik(l,ij)-pbot+rho_bot*(Phi_il(l,ij)-PHI_BOT(ij)) )
+         ml_g2 = gm2*m_il(1,ij)
+         B_il(1,ij) = A_ik(1,ij) + ml_g2 + tau2_g*rho_bot
+         R_il(1,ij) = ml_g2*( Phi_il(1,ij)-Phi_star_il(1,ij)) &
+         + tau2_g*( p_ik(1,ij)-pbot+rho_bot*(Phi_il(1,ij)-PHI_BOT(ij)) )
          DO l = 2,llm
             ! inner interfaces
@@ -68 +67 @@
             ! and residual = tau^2*(ml+(1/g)dl_pi)=0
          END DO
-         l=llm+1
          ! top interface l=llm+1
-         ml_g2 = gm2*m_il(l,ij)
-         B_il(l,ij) = A_ik(l-1,ij) + ml_g2
-         R_il(l,ij) = ml_g2*( Phi_il(l,ij)-Phi_star_il(l,ij)) &
-         + tau2_g*( ptop-p_ik(l-1,ij) )
+         ml_g2 = gm2*m_il(llm+1,ij)
+         B_il(llm+1,ij) = A_ik(llm+1 -1,ij) + ml_g2
+         R_il(llm+1,ij) = ml_g2*( Phi_il(llm+1,ij)-Phi_star_il(llm+1,ij)) &
+         + tau2_g*( ptop-p_ik(llm+1 -1,ij) )
          !
          ! Forward sweep :
          ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
          ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
-         l=1
-         X_ij = 1./B_il(l,ij)
-         C_ik(l,ij) = -A_ik(l,ij) * X_ij
-         D_il(l,ij) = R_il(l,ij) * X_ij
+         X_ij = 1./B_il(1,ij)
+         C_ik(1,ij) = -A_ik(1,ij) * X_ij
+         D_il(1,ij) = R_il(1,ij) * X_ij
          DO l = 2,llm
             X_ij = 1./( B_il(l,ij) + A_ik(l-1,ij)*C_ik(l-1,ij) )
@@ -87 +84 @@
             D_il(l,ij) = (R_il(l,ij)+A_ik(l-1,ij)*D_il(l-1,ij)) * X_ij
          END DO
-         l=llm+1
-         X_ij = 1./( B_il(l,ij) + A_ik(l-1,ij)*C_ik(l-1,ij) )
-         D_il(l,ij) = (R_il(l,ij)+A_ik(l-1,ij)*D_il(l-1,ij)) * X_ij
+         X_ij = 1./( B_il(llm+1,ij) + A_ik(llm+1 -1,ij)*C_ik(llm+1 -1,ij) )
+         D_il(llm+1,ij) = (R_il(llm+1,ij)+A_ik(llm+1 -1,ij)*D_il(llm+1 -1,ij)) * X_ij
          ! Back substitution :
          ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0
          ! + Newton-Raphson update
          ! top interface l=llm+1
-         x_il(l,ij) = D_il(l,ij)
-         Phi_il(l,ij) = Phi_il(l,ij) - x_il(l,ij)
+         x_il(llm+1,ij) = D_il(llm+1,ij)
+         Phi_il(llm+1,ij) = Phi_il(llm+1,ij) - x_il(llm+1,ij)
          DO l = llm,1,-1
             ! Back substitution at lower interfaces