source: codes/icosagcm/devel/src/dynamics/compute_caldyn_slow_NH.F90 @ 876

Last change on this file since 876 was 876, checked in by jisesh, 5 years ago

devel: moved DYSL into compute_caldyn_slow_NH.F90 and compute_caldyn_Coriolis.F90
File size: 9.8 KB
Line 
1MODULE compute_caldyn_slow_NH_mod
2  USE grid_param, ONLY : llm
3  IMPLICIT NONE
4  PRIVATE
5
6#include "../unstructured/unstructured.h90"
7
8  PUBLIC :: compute_caldyn_slow_NH
9
10CONTAINS
11
12#ifdef BEGIN_DYSL
13
14KERNEL(caldyn_slow_NH)
15  FORALL_CELLS_EXT('1', 'llm+1')
16     ON_PRIMAL
17        CST_IF(IS_INNER_INTERFACE, w_il(CELL) = 2.*W(CELL)/(rhodz(KDOWN(CELL))+rhodz(KUP(CELL))) )
18        CST_IF(IS_BOTTOM_LEVEL,    w_il(CELL) = 2.*W(CELL)/rhodz(KUP(CELL)) )
19        CST_IF(IS_TOP_INTERFACE,   w_il(CELL) = 2.*W(CELL)/rhodz(KDOWN(CELL)) )
20     END_BLOCK
21  END_BLOCK
22  FORALL_CELLS_EXT('1', 'llm+1')
23     ON_EDGES
24       ! compute DePhi, v_el, G_el, F_el
25       ! v_el, W2_el and therefore G_el incorporate metric factor le_de
26       ! while DePhil, W_el and F_el do not
27       W_el         = .5*( W(CELL2)+W(CELL1) )
28       DePhil(EDGE) = SIGN*(Phi(CELL2)-Phi(CELL1))
29       F_el(EDGE)   = DePhil(EDGE)*W_el
30       W2_el        = .5*LE_DE * ( W(CELL1)*w_il(CELL1) + W(CELL2)*w_il(CELL2) )
31       v_el(EDGE)   = .5*LE_DE*(u(KUP(EDGE))+u(KDOWN(EDGE))) ! checked
32       G_el(EDGE)   = v_el(EDGE)*W_el - DePhil(EDGE)*W2_el
33     END_BLOCK
34  END_BLOCK
35
36  FORALL_CELLS_EXT('1', 'llm+1')
37    ! compute GradPhi2, dPhi, dW
38    ON_PRIMAL
39       gPhi2=0.
40       dP=0.
41       divG=0
42       FORALL_EDGES   
43         gPhi2 = gPhi2 + LE_DE*DePhil(EDGE)**2
44         dP = dP + LE_DE*DePhil(EDGE)*v_el(EDGE)
45         divG = divG + SIGN*G_el(EDGE) ! -div(G_el), G_el already has le_de
46       END_BLOCK
47       gradPhi2(CELL) = 1./(2.*AI) * gPhi2
48       dPhi(CELL) = gradPhi2(CELL)*w_il(CELL) - 1./(2.*AI)*dP
49       dW(CELL) = (-1./AI)*divG
50    END_BLOCK
51  END_BLOCK
52
53  ! We need a barrier here because we compute gradPhi2, F_el and w_il above and do a vertical average below
54  BARRIER
55
56  FORALL_CELLS_EXT()
57    ! Compute berni at scalar points
58    ON_PRIMAL
59       u2=0.
60       FORALL_EDGES   
61         u2 = u2 + LE_DE*u(EDGE)**2
62       END_BLOCK
63       berni(CELL) = 1./(4.*AI) * u2 - .25*( gradPhi2(CELL)*w_il(CELL)**2 + gradPhi2(UP(CELL))*w_il(UP(CELL))**2 )
64    END_BLOCK
65  END_BLOCK
66
67  FORALL_CELLS_EXT()
68    ON_EDGES
69    ! Compute mass flux and grad(berni)
70      uu = .5*(rhodz(CELL1)+rhodz(CELL2))*u(EDGE) - .5*( F_el(EDGE)+F_el(UP(EDGE)) )
71      hflux(EDGE) = LE_DE*uu
72      du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2))
73    END_BLOCK
74  END_BLOCK
75
76END_BLOCK
77
78#endif END_DYSL
79
80  SUBROUTINE compute_caldyn_vert_NH_unst(mass,geopot,W,wflux, eta_dot,wcov,W_etadot, du,dPhi,dW)
81    USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
82    USE data_unstructured_mod, ONLY : left,right,edge_num,primal_num,dual_num,id_vert_NH, &
83        enter_trace, exit_trace
84    FIELD_MASS   :: mass, eta_dot, wcov, W_etadot ! IN, BUF*3
85    FIELD_GEOPOT :: geopot,W,wflux,dPhi,dW ! IN*3, INOUT*2
86    FIELD_U      :: du ! INOUT
87    DECLARE_INDICES
88    NUM :: w_ij, wflux_ij
89    START_TRACE(id_vert_NH, 6,0,1)
90#include "../kernels_unst/caldyn_vert_NH.k90"
91    STOP_TRACE
92  END SUBROUTINE compute_caldyn_vert_NH_unst
93
94  SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
95    USE icosa
96    USE trace
97    USE caldyn_vars_mod
98    USE omp_para, ONLY : ll_begin, ll_end,ll_beginp1,ll_endp1
99    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
100    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)  ! rho*dz
101    REAL(rstd),INTENT(IN)  :: Phi(iim*jjm,llm+1)  ! prognostic geopotential
102    REAL(rstd),INTENT(IN)  :: W(iim*jjm,llm+1)    ! prognostic vertical momentum
103
104    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
105    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)
106    REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1)
107    REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1)
108
109    REAL(rstd) :: w_il(iim*jjm,llm+1) ! Wil/mil
110    REAL(rstd) :: F_el(3*iim*jjm,llm+1) ! NH mass flux
111    REAL(rstd) :: gradPhi2(iim*jjm,llm+1) ! grad_Phi**2
112    REAL(rstd) :: DePhil(3*iim*jjm,llm+1) ! grad(Phi)
113   
114    INTEGER :: ij,l,kdown,kup
115    REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu
116
117    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
118    REAL(rstd) :: G_el(3*iim*jjm,llm+1) ! horizontal flux of W
119    REAL(rstd) :: v_el(3*iim*jjm,llm+1)
120
121    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
122    REAL(rstd) :: G_el1(3*iim*jjm) ! horizontal flux of W
123    REAL(rstd) :: v_el1(3*iim*jjm)
124
125    CALL trace_start("compute_caldyn_slow_NH")
126
127    IF(dysl) THEN
128
129!$OMP BARRIER
130#include "../kernels_hex/caldyn_slow_NH.k90"
131!$OMP BARRIER
132     
133     ELSE
134
135#define BERNI(ij) berni1(ij)
136#define G_EL(ij) G_el1(ij)
137#define V_EL(ij) v_el1(ij)
138
139    DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
140       IF(l==1) THEN
141          kdown=1
142       ELSE
143          kdown=l-1
144       END IF
145       IF(l==llm+1) THEN
146          kup=llm
147       ELSE
148          kup=l
149       END IF
150       ! below : "checked" means "formula also valid when kup=kdown (top/bottom)"
151       ! compute mil, wil=Wil/mil
152       DO ij=ij_begin_ext, ij_end_ext
153          w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) ! checked
154       END DO
155       ! compute DePhi, v_el, G_el, F_el
156       ! v_el, W2_el and therefore G_el incorporate metric factor le_de
157       ! while DePhil, W_el and F_el don't
158       DO ij=ij_begin_ext, ij_end_ext
159          ! Compute on edge 'right'
160          W_el  = .5*( W(ij,l)+W(ij+t_right,l) )
161          DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
162          F_el(ij+u_right,l)   = DePhil(ij+u_right,l)*W_el
163          W2_el = .5*le_de(ij+u_right) * &
164               ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) )
165          V_EL(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked
166          G_EL(ij+u_right) = V_EL(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el
167          ! Compute on edge 'lup'
168          W_el  = .5*( W(ij,l)+W(ij+t_lup,l) )
169          DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
170          F_el(ij+u_lup,l)   = DePhil(ij+u_lup,l)*W_el
171          W2_el = .5*le_de(ij+u_lup) * &
172               ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) )
173          V_EL(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked
174          G_EL(ij+u_lup) = V_EL(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el
175          ! Compute on edge 'ldown'
176          W_el  = .5*( W(ij,l)+W(ij+t_ldown,l) )
177          DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
178          F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el
179          W2_el = .5*le_de(ij+u_ldown) * &
180               ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) )
181          V_EL(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked
182          G_EL(ij+u_ldown) = V_EL(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el
183       END DO
184       ! compute GradPhi2, dPhi, dW
185       DO ij=ij_begin_ext, ij_end_ext
186          gradPhi2(ij,l) = &
187               1/(2*Ai(ij))*(le_de(ij+u_right)*DePhil(ij+u_right,l)**2 + &
188               le_de(ij+u_rup)*DePhil(ij+u_rup,l)**2 +      &
189               le_de(ij+u_lup)*DePhil(ij+u_lup,l)**2 +      &
190               le_de(ij+u_left)*DePhil(ij+u_left,l)**2 +    &
191               le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 +  &
192               le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 )
193
194          dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* & 
195               ( DePhil(ij+u_right,l)*V_EL(ij+u_right) + & ! -v.gradPhi,
196                 DePhil(ij+u_rup,l)*V_EL(ij+u_rup) +     & ! v_el already has le_de
197                 DePhil(ij+u_lup,l)*V_EL(ij+u_lup) +     &
198                 DePhil(ij+u_left,l)*V_EL(ij+u_left) +   &
199                 DePhil(ij+u_ldown,l)*V_EL(ij+u_ldown) + &
200                 DePhil(ij+u_rdown,l)*V_EL(ij+u_rdown) )
201
202          dW(ij,l) = -1./Ai(ij)*(           & ! -div(G_el),
203               ne_right*G_EL(ij+u_right) +  & ! G_el already has le_de
204               ne_rup*G_EL(ij+u_rup) +      &
205               ne_lup*G_EL(ij+u_lup) +      & 
206               ne_left*G_EL(ij+u_left) +    &
207               ne_ldown*G_EL(ij+u_ldown) +  &
208               ne_rdown*G_EL(ij+u_rdown))
209       END DO
210    END DO
211
212    DO l=ll_begin, ll_end ! compute on k levels (layers)
213       ! Compute berni at scalar points
214       DO ij=ij_begin_ext, ij_end_ext
215          BERNI(ij) = &
216               1/(4*Ai(ij))*( &
217                    le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
218                    le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
219                    le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
220                    le_de(ij+u_left)*u(ij+u_left,l)**2 +       &
221                    le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
222                    le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) &
223               - .25*( gradPhi2(ij,l)  *w_il(ij,l)**2 +   &
224                      gradPhi2(ij,l+1)*w_il(ij,l+1)**2 )
225       END DO
226       ! Compute mass flux and grad(berni) at edges
227       DO ij=ij_begin_ext, ij_end_ext
228          ! Compute on edge 'right'
229          uu_right = 0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) &
230                    -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1))
231          hflux(ij+u_right,l) = uu_right*le_de(ij+u_right)
232          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
233          ! Compute on edge 'lup'
234          uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) &
235                  -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))
236          hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup)
237          du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
238          ! Compute on edge 'ldown'
239          uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) &
240                    -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1))
241          hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown)
242          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
243       END DO
244    END DO
245
246#undef V_EL
247#undef G_EL
248#undef BERNI
249
250    END IF ! dysl
251
252    CALL trace_end("compute_caldyn_slow_NH")
253
254  END SUBROUTINE compute_caldyn_slow_NH
255
256END MODULE compute_caldyn_slow_NH_mod
Note: See TracBrowser for help on using the repository browser.