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

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

devel: included USE module trace in compute_caldyn_slow_NH and compute_caldyn_slow_hydro

File size: 7.0 KB
Line 
1MODULE compute_caldyn_slow_NH_mod
2  USE grid_param, ONLY : llm
3  IMPLICIT NONE
4  PRIVATE
5
6  PUBLIC :: compute_caldyn_slow_NH
7
8CONTAINS
9
10  SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
11    USE icosa
12    USE trace
13    USE caldyn_vars_mod
14    USE omp_para, ONLY : ll_begin, ll_end,ll_beginp1,ll_endp1
15    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
16    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)  ! rho*dz
17    REAL(rstd),INTENT(IN)  :: Phi(iim*jjm,llm+1)  ! prognostic geopotential
18    REAL(rstd),INTENT(IN)  :: W(iim*jjm,llm+1)    ! prognostic vertical momentum
19
20    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
21    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)
22    REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1)
23    REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1)
24
25    REAL(rstd) :: w_il(iim*jjm,llm+1) ! Wil/mil
26    REAL(rstd) :: F_el(3*iim*jjm,llm+1) ! NH mass flux
27    REAL(rstd) :: gradPhi2(iim*jjm,llm+1) ! grad_Phi**2
28    REAL(rstd) :: DePhil(3*iim*jjm,llm+1) ! grad(Phi)
29   
30    INTEGER :: ij,l,kdown,kup
31    REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu
32
33    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
34    REAL(rstd) :: G_el(3*iim*jjm,llm+1) ! horizontal flux of W
35    REAL(rstd) :: v_el(3*iim*jjm,llm+1)
36
37    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
38    REAL(rstd) :: G_el1(3*iim*jjm) ! horizontal flux of W
39    REAL(rstd) :: v_el1(3*iim*jjm)
40
41    CALL trace_start("compute_caldyn_slow_NH")
42
43    IF(dysl) THEN
44
45!$OMP BARRIER
46#include "../kernels_hex/caldyn_slow_NH.k90"
47!$OMP BARRIER
48     
49     ELSE
50
51#define BERNI(ij) berni1(ij)
52#define G_EL(ij) G_el1(ij)
53#define V_EL(ij) v_el1(ij)
54
55    DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
56       IF(l==1) THEN
57          kdown=1
58       ELSE
59          kdown=l-1
60       END IF
61       IF(l==llm+1) THEN
62          kup=llm
63       ELSE
64          kup=l
65       END IF
66       ! below : "checked" means "formula also valid when kup=kdown (top/bottom)"
67       ! compute mil, wil=Wil/mil
68       DO ij=ij_begin_ext, ij_end_ext
69          w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) ! checked
70       END DO
71       ! compute DePhi, v_el, G_el, F_el
72       ! v_el, W2_el and therefore G_el incorporate metric factor le_de
73       ! while DePhil, W_el and F_el don't
74       DO ij=ij_begin_ext, ij_end_ext
75          ! Compute on edge 'right'
76          W_el  = .5*( W(ij,l)+W(ij+t_right,l) )
77          DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
78          F_el(ij+u_right,l)   = DePhil(ij+u_right,l)*W_el
79          W2_el = .5*le_de(ij+u_right) * &
80               ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) )
81          V_EL(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked
82          G_EL(ij+u_right) = V_EL(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el
83          ! Compute on edge 'lup'
84          W_el  = .5*( W(ij,l)+W(ij+t_lup,l) )
85          DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
86          F_el(ij+u_lup,l)   = DePhil(ij+u_lup,l)*W_el
87          W2_el = .5*le_de(ij+u_lup) * &
88               ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) )
89          V_EL(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked
90          G_EL(ij+u_lup) = V_EL(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el
91          ! Compute on edge 'ldown'
92          W_el  = .5*( W(ij,l)+W(ij+t_ldown,l) )
93          DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
94          F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el
95          W2_el = .5*le_de(ij+u_ldown) * &
96               ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) )
97          V_EL(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked
98          G_EL(ij+u_ldown) = V_EL(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el
99       END DO
100       ! compute GradPhi2, dPhi, dW
101       DO ij=ij_begin_ext, ij_end_ext
102          gradPhi2(ij,l) = &
103               1/(2*Ai(ij))*(le_de(ij+u_right)*DePhil(ij+u_right,l)**2 + &
104               le_de(ij+u_rup)*DePhil(ij+u_rup,l)**2 +      &
105               le_de(ij+u_lup)*DePhil(ij+u_lup,l)**2 +      &
106               le_de(ij+u_left)*DePhil(ij+u_left,l)**2 +    &
107               le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 +  &
108               le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 )
109
110          dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* & 
111               ( DePhil(ij+u_right,l)*V_EL(ij+u_right) + & ! -v.gradPhi,
112                 DePhil(ij+u_rup,l)*V_EL(ij+u_rup) +     & ! v_el already has le_de
113                 DePhil(ij+u_lup,l)*V_EL(ij+u_lup) +     &
114                 DePhil(ij+u_left,l)*V_EL(ij+u_left) +   &
115                 DePhil(ij+u_ldown,l)*V_EL(ij+u_ldown) + &
116                 DePhil(ij+u_rdown,l)*V_EL(ij+u_rdown) )
117
118          dW(ij,l) = -1./Ai(ij)*(           & ! -div(G_el),
119               ne_right*G_EL(ij+u_right) +  & ! G_el already has le_de
120               ne_rup*G_EL(ij+u_rup) +      &
121               ne_lup*G_EL(ij+u_lup) +      & 
122               ne_left*G_EL(ij+u_left) +    &
123               ne_ldown*G_EL(ij+u_ldown) +  &
124               ne_rdown*G_EL(ij+u_rdown))
125       END DO
126    END DO
127
128    DO l=ll_begin, ll_end ! compute on k levels (layers)
129       ! Compute berni at scalar points
130       DO ij=ij_begin_ext, ij_end_ext
131          BERNI(ij) = &
132               1/(4*Ai(ij))*( &
133                    le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
134                    le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
135                    le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
136                    le_de(ij+u_left)*u(ij+u_left,l)**2 +       &
137                    le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
138                    le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) &
139               - .25*( gradPhi2(ij,l)  *w_il(ij,l)**2 +   &
140                      gradPhi2(ij,l+1)*w_il(ij,l+1)**2 )
141       END DO
142       ! Compute mass flux and grad(berni) at edges
143       DO ij=ij_begin_ext, ij_end_ext
144          ! Compute on edge 'right'
145          uu_right = 0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) &
146                    -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1))
147          hflux(ij+u_right,l) = uu_right*le_de(ij+u_right)
148          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
149          ! Compute on edge 'lup'
150          uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) &
151                  -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))
152          hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup)
153          du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
154          ! Compute on edge 'ldown'
155          uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) &
156                    -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1))
157          hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown)
158          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
159       END DO
160    END DO
161
162#undef V_EL
163#undef G_EL
164#undef BERNI
165
166    END IF ! dysl
167
168    CALL trace_end("compute_caldyn_slow_NH")
169
170  END SUBROUTINE compute_caldyn_slow_NH
171
172END MODULE compute_caldyn_slow_NH_mod
Note: See TracBrowser for help on using the repository browser.