source: codes/icosagcm/devel/DySL/kernels_caldyn_NH.jin @ 941

Last change on this file since 941 was 941, checked in by dubos, 5 years ago

devel : split DySL from Python
File size: 2.0 KB
Line 
1{% macro compute_p_and_Aik() %}
2{% set p_and_c2_from_theta=caller %}
3    SEQUENCE_C1
4      BODY('1,llm')
5          rho_ij = (g*m_ik(CELL))/(Phi_il(UP(CELL))-Phi_il(CELL))
6          {{ p_and_c2_from_theta() }}
7          A_ik(CELL) = c2_mik*(tau/g*rho_ij)**2
8      END_BLOCK
9    END_BLOCK
10{%- endmacro %}
11
12KERNEL(caldyn_mil)
13  FORALL_CELLS_EXT('1', 'llm+1')
14    ON_PRIMAL
15       CST_IF(IS_BOTTOM_LEVEL,    m_il(CELL) = .5*rhodz(CELL)                     )
16       CST_IF(IS_INNER_INTERFACE, m_il(CELL) = .5*(rhodz(CELL)+rhodz(DOWN(CELL))) )
17       CST_IF(IS_TOP_INTERFACE,   m_il(CELL) = .5*rhodz(DOWN(CELL))               )
18    END_BLOCK
19  END_BLOCK
20END_BLOCK
21
22KERNEL(caldyn_vert_NH)
23  FORALL_CELLS_EXT()
24    ON_PRIMAL
25      CST_IF(IS_BOTTOM_LEVEL,     w_ij = ( W(CELL)+.5*W(UP(CELL)) )/mass(CELL) )
26      CST_IF(IS_INNER_LAYER , w_ij = .5*( W(CELL)+W(UP(CELL)) )/mass(CELL) )
27      CST_IF(IS_TOP_LAYER,        w_ij = ( .5*W(CELL)+W(UP(CELL)) )/mass(CELL) )
28      wflux_ij = .5*(wflux(CELL)+wflux(UP(CELL)))
29      W_etadot(CELL) = wflux_ij*w_ij
30      eta_dot(CELL) = wflux_ij / mass(CELL)
31      wcov(CELL) = w_ij*(geopot(UP(CELL))-geopot(CELL))
32    END_BLOCK
33  END_BLOCK
34  ! add NH term to du
35  FORALL_CELLS()
36    ON_EDGES
37      du(EDGE) = du(EDGE) - .5*(wcov(CELL1)+wcov(CELL2))*SIGN*(eta_dot(CELL2)-eta_dot(CELL1))
38    END_BLOCK
39  END_BLOCK
40  ! add NH terms to dW, dPhi
41  ! FIXME : TODO top and bottom
42  FORALL_CELLS('2','llm')
43    ON_PRIMAL
44      dPhi(CELL)=dPhi(CELL)-wflux(CELL)*(geopot(UP(CELL))-geopot(DOWN(CELL)))/(mass(DOWN(CELL))+mass(CELL))
45    END_BLOCK
46  END_BLOCK
47
48  ! We need a barrier here because we compute W_etadot above and do a vertical difference below
49  BARRIER
50
51  FORALL_CELLS('1','llm+1')
52    ON_PRIMAL
53      CST_IFTHEN(IS_BOTTOM_LEVEL)
54         dW(CELL) = dW(CELL) - W_etadot(CELL)
55      CST_ELSEIF(IS_TOP_INTERFACE)
56         dW(CELL) = dW(CELL) + W_etadot(DOWN(CELL))
57      CST_ELSE
58         dW(CELL) = dW(CELL) + W_etadot(DOWN(CELL)) - W_etadot(CELL)
59      CST_ENDIF
60    END_BLOCK
61  END_BLOCK
62END_BLOCK
63
Note: See TracBrowser for help on using the repository browser.