Changeset 686 for codes/icosagcm

Timestamp:

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

Author:

dubos

Message:

devel/unstructured : piecewise-constant vertical remapping

Location:

codes/icosagcm/devel/src

Files:

• icosa_init.f90 (modified) (1 diff)
• kernels_unst/caldyn_wflux.k90 (modified) (1 diff)
• kernels_unst/compute_NH_geopot.k90 (modified) (3 diffs)
• kernels_unst/compute_geopot.k90 (modified) (4 diffs)
• kernels_unst/remap_eta.k90 (added)
• kernels_unst/remap_theta.k90 (added)
• kernels_unst/remap_u.k90 (added)
• kernels_unst/theta.k90 (modified) (1 diff)
• unstructured/transport_unstructured.F90 (added)

codes/icosagcm/devel/src/icosa_init.f90

@@ -109 +109 @@
   SUBROUTINE force_recompile_unstructured
     USE timestep_unstructured_mod
+    USE transport_unstructured_mod
   END SUBROUTINE force_recompile_unstructured
 

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

@@ -7 +7 @@
          convm(l,ij) = convm(l,ij) + convm(l+1,ij)
       END DO
-      l=1
-      dmass_col(ij) = convm(l,ij)
+      dmass_col(ij) = convm(1,ij)
       DO l = 2,llm
          ! Compute vertical mass flux (l=1,llm+1 set to zero at init)

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

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

@@ -6 +6 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij=1,primal_num
-         l=llm
-         pk(l,ij) = ptop + .5*g* theta(l,ij,1)*rhodz(l,ij)
+         pk(llm,ij) = ptop + .5*g* theta(llm,ij,1)*rhodz(llm,ij)
          DO l = llm-1,1,-1
             pk(l,ij) = pk(l+1,ij) + (.5*g)*( theta(l,ij,1)*rhodz(l,ij) + theta(l+1,ij,1)*rhodz(l+1,ij) )
          END DO
          IF(caldyn_eta == eta_lag) THEN
-            l=1
-            ps(ij) = pk(l,ij) + .5*g* theta(l,ij,1)*rhodz(l,ij)
+            ps(ij) = pk(1,ij) + .5*g* theta(1,ij,1)*rhodz(1,ij)
          END IF
       END DO
@@ -29 +27 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij=1,primal_num
-         l=llm
-         pk(l,ij) = ptop + .5*g* rhodz(l,ij)
+         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)
          DO l = llm-1,1,-1
             pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij) + rhodz(l+1,ij) )
          END DO
          IF(caldyn_eta == eta_lag) THEN
-            l=1
-            ps(ij) = pk(l,ij) + .5*g* rhodz(l,ij)
+            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)
          END IF
       END DO
@@ -54 +50 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij=1,primal_num
-         l=llm
-         pk(l,ij) = ptop + .5*g* rhodz(l,ij)
+         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)
          DO l = llm-1,1,-1
             pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij) + rhodz(l+1,ij) )
          END DO
          IF(caldyn_eta == eta_lag) THEN
-            l=1
-            ps(ij) = pk(l,ij) + .5*g* rhodz(l,ij)
+            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)
          END IF
       END DO
@@ -79 +73 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij=1,primal_num
-         l=llm
-         pk(l,ij) = ptop + .5*g* rhodz(l,ij)*(1.+theta(l,ij,2))
+         pk(llm,ij) = ptop + .5*g* rhodz(llm,ij)*(1.+theta(llm,ij,2))
          DO l = llm-1,1,-1
             pk(l,ij) = pk(l+1,ij) + (.5*g)*( rhodz(l,ij)*(1.+theta(l,ij,2)) + rhodz(l+1,ij)*(1.+theta(l+1,ij,2)) )
          END DO
          IF(caldyn_eta == eta_lag) THEN
-            l=1
-            ps(ij) = pk(l,ij) + .5*g* rhodz(l,ij)*(1.+theta(l,ij,2))
+            ps(ij) = pk(1,ij) + .5*g* rhodz(1,ij)*(1.+theta(1,ij,2))
          END IF
       END DO

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

@@ -7 +7 @@
       !$OMP DO SCHEDULE(STATIC)
       DO ij=1,primal_num
-         l=0
          mass_col(ij)=0.
          DO l = 1,llm