Changeset 876 for codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_NH.F90

Timestamp:

05/28/19 22:44:58 (5 years ago)

Author:

jisesh

Message:

devel: moved DYSL into compute_caldyn_slow_NH.F90 and compute_caldyn_Coriolis.F90

File:

• codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_NH.F90 (modified) (1 diff)

codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_NH.F90

@@ -4 +4 @@
   PRIVATE
 
+#include "../unstructured/unstructured.h90"
+
   PUBLIC :: compute_caldyn_slow_NH
 
 CONTAINS
+
+#ifdef BEGIN_DYSL
+
+KERNEL(caldyn_slow_NH)
+  FORALL_CELLS_EXT('1', 'llm+1')
+     ON_PRIMAL
+        CST_IF(IS_INNER_INTERFACE, w_il(CELL) = 2.*W(CELL)/(rhodz(KDOWN(CELL))+rhodz(KUP(CELL))) )
+        CST_IF(IS_BOTTOM_LEVEL,    w_il(CELL) = 2.*W(CELL)/rhodz(KUP(CELL)) )
+        CST_IF(IS_TOP_INTERFACE,   w_il(CELL) = 2.*W(CELL)/rhodz(KDOWN(CELL)) )
+     END_BLOCK
+  END_BLOCK
+  FORALL_CELLS_EXT('1', 'llm+1')
+     ON_EDGES
+       ! compute DePhi, v_el, G_el, F_el
+       ! v_el, W2_el and therefore G_el incorporate metric factor le_de
+       ! while DePhil, W_el and F_el do not 
+       W_el         = .5*( W(CELL2)+W(CELL1) )
+       DePhil(EDGE) = SIGN*(Phi(CELL2)-Phi(CELL1))
+       F_el(EDGE)   = DePhil(EDGE)*W_el
+       W2_el        = .5*LE_DE * ( W(CELL1)*w_il(CELL1) + W(CELL2)*w_il(CELL2) )
+       v_el(EDGE)   = .5*LE_DE*(u(KUP(EDGE))+u(KDOWN(EDGE))) ! checked
+       G_el(EDGE)   = v_el(EDGE)*W_el - DePhil(EDGE)*W2_el
+     END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS_EXT('1', 'llm+1')
+    ! compute GradPhi2, dPhi, dW
+    ON_PRIMAL
+       gPhi2=0.
+       dP=0.
+       divG=0
+       FORALL_EDGES    
+         gPhi2 = gPhi2 + LE_DE*DePhil(EDGE)**2
+         dP = dP + LE_DE*DePhil(EDGE)*v_el(EDGE)
+         divG = divG + SIGN*G_el(EDGE) ! -div(G_el), G_el already has le_de
+       END_BLOCK
+       gradPhi2(CELL) = 1./(2.*AI) * gPhi2
+       dPhi(CELL) = gradPhi2(CELL)*w_il(CELL) - 1./(2.*AI)*dP
+       dW(CELL) = (-1./AI)*divG
+    END_BLOCK
+  END_BLOCK
+
+  ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below
+  BARRIER
+
+  FORALL_CELLS_EXT()
+    ! Compute berni at scalar points
+    ON_PRIMAL
+       u2=0.
+       FORALL_EDGES    
+         u2 = u2 + LE_DE*u(EDGE)**2
+       END_BLOCK
+       berni(CELL) = 1./(4.*AI) * u2 - .25*( gradPhi2(CELL)*w_il(CELL)**2 + gradPhi2(UP(CELL))*w_il(UP(CELL))**2 )
+    END_BLOCK
+  END_BLOCK
+
+  FORALL_CELLS_EXT()
+    ON_EDGES
+    ! Compute mass flux and grad(berni)
+      uu = .5*(rhodz(CELL1)+rhodz(CELL2))*u(EDGE) - .5*( F_el(EDGE)+F_el(UP(EDGE)) )
+      hflux(EDGE) = LE_DE*uu
+      du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2))
+    END_BLOCK
+  END_BLOCK
+
+END_BLOCK
+
+#endif END_DYSL
+
+  SUBROUTINE compute_caldyn_vert_NH_unst(mass,geopot,W,wflux, eta_dot,wcov,W_etadot, du,dPhi,dW)
+    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
+    USE data_unstructured_mod, ONLY : left,right,edge_num,primal_num,dual_num,id_vert_NH, &
+        enter_trace, exit_trace
+    FIELD_MASS   :: mass, eta_dot, wcov, W_etadot ! IN, BUF*3
+    FIELD_GEOPOT :: geopot,W,wflux,dPhi,dW ! IN*3, INOUT*2
+    FIELD_U      :: du ! INOUT
+    DECLARE_INDICES
+    NUM :: w_ij, wflux_ij
+    START_TRACE(id_vert_NH, 6,0,1)
+#include "../kernels_unst/caldyn_vert_NH.k90"
+    STOP_TRACE
+  END SUBROUTINE compute_caldyn_vert_NH_unst
 
   SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)