Index: codes/icosagcm/devel/src/kernels_unst/wind_centered.k90
===================================================================
--- codes/icosagcm/devel/src/kernels_unst/wind_centered.k90	(revision 796)
+++ codes/icosagcm/devel/src/kernels_unst/wind_centered.k90	(revision 796)
@@ -0,0 +1,190 @@
+   !--------------------------------------------------------------------------
+   !---------------------------- wind_centered ----------------------------------
+   ! Perot reconstruction based on Gauss theorem
+   ! u = sum( u.edge_normal * edge_length * (edge_midpoint-cell_centroid) ) /cell_area
+   !$OMP DO SCHEDULE(STATIC)
+   DO ij = 1, primal_num
+      ! this VLOOP iterates over primal cell edges
+      SELECT CASE(primal_deg(ij))
+      CASE(4)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         ij_up1 = up(edge1)
+         ij_up2 = up(edge2)
+         ij_up3 = up(edge3)
+         ij_up4 = up(edge4)
+         ij_down1 = down(edge1)
+         ij_down2 = down(edge2)
+         ij_down3 = down(edge3)
+         ij_down4 = down(edge4)
+         !DIR$ SIMD
+         DO l = 1, llm
+            cx=centroid(ij,1)
+            cy=centroid(ij,2)
+            cz=centroid(ij,3)
+            ux=0. ; uy=0. ; uz=0.
+            ue_le = sign1*ue(l,edge1)*le(edge1)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up1,1)+xyz_v(ij_down1,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up1,2)+xyz_v(ij_down1,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up1,3)+xyz_v(ij_down1,3))-cz)
+            ue_le = sign2*ue(l,edge2)*le(edge2)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up2,1)+xyz_v(ij_down2,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up2,2)+xyz_v(ij_down2,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up2,3)+xyz_v(ij_down2,3))-cz)
+            ue_le = sign3*ue(l,edge3)*le(edge3)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up3,1)+xyz_v(ij_down3,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up3,2)+xyz_v(ij_down3,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up3,3)+xyz_v(ij_down3,3))-cz)
+            ue_le = sign4*ue(l,edge4)*le(edge4)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up4,1)+xyz_v(ij_down4,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up4,2)+xyz_v(ij_down4,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up4,3)+xyz_v(ij_down4,3))-cz)
+            fac = scale*(1./Ai(ij))
+            ucenter(l,ij,1)=ux*fac
+            ucenter(l,ij,2)=uy*fac
+            ucenter(l,ij,3)=uz*fac
+         END DO
+      CASE(5)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         sign5 = primal_ne(5,ij)
+         ij_up1 = up(edge1)
+         ij_up2 = up(edge2)
+         ij_up3 = up(edge3)
+         ij_up4 = up(edge4)
+         ij_up5 = up(edge5)
+         ij_down1 = down(edge1)
+         ij_down2 = down(edge2)
+         ij_down3 = down(edge3)
+         ij_down4 = down(edge4)
+         ij_down5 = down(edge5)
+         !DIR$ SIMD
+         DO l = 1, llm
+            cx=centroid(ij,1)
+            cy=centroid(ij,2)
+            cz=centroid(ij,3)
+            ux=0. ; uy=0. ; uz=0.
+            ue_le = sign1*ue(l,edge1)*le(edge1)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up1,1)+xyz_v(ij_down1,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up1,2)+xyz_v(ij_down1,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up1,3)+xyz_v(ij_down1,3))-cz)
+            ue_le = sign2*ue(l,edge2)*le(edge2)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up2,1)+xyz_v(ij_down2,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up2,2)+xyz_v(ij_down2,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up2,3)+xyz_v(ij_down2,3))-cz)
+            ue_le = sign3*ue(l,edge3)*le(edge3)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up3,1)+xyz_v(ij_down3,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up3,2)+xyz_v(ij_down3,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up3,3)+xyz_v(ij_down3,3))-cz)
+            ue_le = sign4*ue(l,edge4)*le(edge4)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up4,1)+xyz_v(ij_down4,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up4,2)+xyz_v(ij_down4,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up4,3)+xyz_v(ij_down4,3))-cz)
+            ue_le = sign5*ue(l,edge5)*le(edge5)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up5,1)+xyz_v(ij_down5,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up5,2)+xyz_v(ij_down5,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up5,3)+xyz_v(ij_down5,3))-cz)
+            fac = scale*(1./Ai(ij))
+            ucenter(l,ij,1)=ux*fac
+            ucenter(l,ij,2)=uy*fac
+            ucenter(l,ij,3)=uz*fac
+         END DO
+      CASE(6)
+         edge1 = primal_edge(1,ij)
+         edge2 = primal_edge(2,ij)
+         edge3 = primal_edge(3,ij)
+         edge4 = primal_edge(4,ij)
+         edge5 = primal_edge(5,ij)
+         edge6 = primal_edge(6,ij)
+         sign1 = primal_ne(1,ij)
+         sign2 = primal_ne(2,ij)
+         sign3 = primal_ne(3,ij)
+         sign4 = primal_ne(4,ij)
+         sign5 = primal_ne(5,ij)
+         sign6 = primal_ne(6,ij)
+         ij_up1 = up(edge1)
+         ij_up2 = up(edge2)
+         ij_up3 = up(edge3)
+         ij_up4 = up(edge4)
+         ij_up5 = up(edge5)
+         ij_up6 = up(edge6)
+         ij_down1 = down(edge1)
+         ij_down2 = down(edge2)
+         ij_down3 = down(edge3)
+         ij_down4 = down(edge4)
+         ij_down5 = down(edge5)
+         ij_down6 = down(edge6)
+         !DIR$ SIMD
+         DO l = 1, llm
+            cx=centroid(ij,1)
+            cy=centroid(ij,2)
+            cz=centroid(ij,3)
+            ux=0. ; uy=0. ; uz=0.
+            ue_le = sign1*ue(l,edge1)*le(edge1)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up1,1)+xyz_v(ij_down1,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up1,2)+xyz_v(ij_down1,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up1,3)+xyz_v(ij_down1,3))-cz)
+            ue_le = sign2*ue(l,edge2)*le(edge2)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up2,1)+xyz_v(ij_down2,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up2,2)+xyz_v(ij_down2,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up2,3)+xyz_v(ij_down2,3))-cz)
+            ue_le = sign3*ue(l,edge3)*le(edge3)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up3,1)+xyz_v(ij_down3,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up3,2)+xyz_v(ij_down3,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up3,3)+xyz_v(ij_down3,3))-cz)
+            ue_le = sign4*ue(l,edge4)*le(edge4)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up4,1)+xyz_v(ij_down4,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up4,2)+xyz_v(ij_down4,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up4,3)+xyz_v(ij_down4,3))-cz)
+            ue_le = sign5*ue(l,edge5)*le(edge5)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up5,1)+xyz_v(ij_down5,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up5,2)+xyz_v(ij_down5,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up5,3)+xyz_v(ij_down5,3))-cz)
+            ue_le = sign6*ue(l,edge6)*le(edge6)
+            ux = ux + ue_le*(.5*(xyz_v(ij_up6,1)+xyz_v(ij_down6,1))-cx)
+            uy = uy + ue_le*(.5*(xyz_v(ij_up6,2)+xyz_v(ij_down6,2))-cy)
+            uz = uz + ue_le*(.5*(xyz_v(ij_up6,3)+xyz_v(ij_down6,3))-cz)
+            fac = scale*(1./Ai(ij))
+            ucenter(l,ij,1)=ux*fac
+            ucenter(l,ij,2)=uy*fac
+            ucenter(l,ij,3)=uz*fac
+         END DO
+      CASE DEFAULT
+         !DIR$ SIMD
+         DO l = 1, llm
+            cx=centroid(ij,1)
+            cy=centroid(ij,2)
+            cz=centroid(ij,3)
+            ux=0. ; uy=0. ; uz=0.
+            DO iedge = 1, primal_deg(ij)
+               edge = primal_edge(iedge,ij)
+               ij_up = up(edge)
+               ij_down = down(edge)
+               ue_le = primal_ne(iedge,ij)*ue(l,edge)*le(edge)
+               ux = ux + ue_le*(.5*(xyz_v(ij_up,1)+xyz_v(ij_down,1))-cx)
+               uy = uy + ue_le*(.5*(xyz_v(ij_up,2)+xyz_v(ij_down,2))-cy)
+               uz = uz + ue_le*(.5*(xyz_v(ij_up,3)+xyz_v(ij_down,3))-cz)
+            END DO
+            fac = scale*(1./Ai(ij))
+            ucenter(l,ij,1)=ux*fac
+            ucenter(l,ij,2)=uy*fac
+            ucenter(l,ij,3)=uz*fac
+         END DO
+      END SELECT
+   END DO
+   !$OMP END DO
+   !---------------------------- wind_centered ----------------------------------
+   !--------------------------------------------------------------------------
