source: codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90 @ 842

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

devel: separate module for compute_caldyn_kv
File size: 39.3 KB
Line 
1MODULE caldyn_kernels_hevi_mod
2  USE icosa
3  USE trace
4  USE omp_para
5  USE disvert_mod
6  USE transfert_mod
7  USE caldyn_vars_mod
8  IMPLICIT NONE
9  PRIVATE
10
11  REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6
12
13  LOGICAL, SAVE :: debug_hevi_solver = .FALSE.
14
15  PUBLIC :: compute_caldyn_Coriolis, &
16       compute_caldyn_slow_hydro, compute_caldyn_slow_NH, &
17       compute_caldyn_solver, compute_caldyn_fast
18
19CONTAINS
20
21  SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il)
22    REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs
23    REAL(rstd),INTENT(IN)    :: phis(iim*jjm)
24    REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm)
25    REAL(rstd),INTENT(IN)    :: m_il(iim*jjm,llm+1)
26    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
27    REAL(rstd),INTENT(IN)    :: W_il(iim*jjm,llm+1)   ! vertical momentum
28    REAL(rstd),INTENT(INOUT) :: Phi_il(iim*jjm,llm+1) ! geopotential
29
30    REAL(rstd) :: Phi_star_il(iim*jjm,llm+1) 
31    REAL(rstd) :: p_ik(iim*jjm,llm)      ! pressure
32    REAL(rstd) :: R_il(iim*jjm,llm+1)    ! rhs of tridiag problem
33    REAL(rstd) :: x_il(iim*jjm,llm+1)    ! solution of tridiag problem
34    REAL(rstd) :: A_ik(iim*jjm,llm)      ! off-diagonal coefficients of tridiag problem
35    REAL(rstd) :: B_il(iim*jjm,llm+1)    ! diagonal coefficients of tridiag problem
36    REAL(rstd) :: C_ik(iim*jjm,llm)      ! Thomas algorithm
37    REAL(rstd) :: D_il(iim*jjm,llm+1)    ! Thomas algorithm
38    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij
39    REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik, vreff
40
41    INTEGER    :: iter, ij, l, ij_omp_begin_ext, ij_omp_end_ext
42
43    CALL distrib_level(ij_begin_ext,ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext)
44
45    IF(dysl) THEN
46#define PHI_BOT(ij) phis(ij)
47#include "../kernels_hex/compute_NH_geopot.k90"
48#undef PHI_BOT
49    ELSE
50!    FIXME : vertical OpenMP parallelism will not work
51   
52    tau2_g=tau*tau/g
53    g2=g*g
54    gm2 = g**-2
55    gamma = 1./(1.-kappa)
56   
57    ! compute Phi_star
58    DO l=1,llm+1
59       !DIR$ SIMD
60       DO ij=ij_begin_ext,ij_end_ext
61          Phi_star_il(ij,l) = Phi_il(ij,l) + tau*g2*(W_il(ij,l)/m_il(ij,l)-tau)
62       ENDDO
63    ENDDO
64
65    ! Newton-Raphson iteration : Phi_il contains current guess value
66    DO iter=1,5 ! 2 iterations should be enough
67
68       ! Compute pressure, A_ik
69       DO l=1,llm
70          !DIR$ SIMD
71          DO ij=ij_begin_ext,ij_end_ext
72             rho_ij = (g*m_ik(ij,l))/(Phi_il(ij,l+1)-Phi_il(ij,l))
73             X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij
74             p_ik(ij,l) = preff*(X_ij**gamma)
75             c2_mik = gamma*p_ik(ij,l)/(rho_ij*m_ik(ij,l)) ! c^2 = gamma*R*T = gamma*p/rho
76             A_ik(ij,l) = c2_mik*(tau/g*rho_ij)**2
77          ENDDO
78       ENDDO
79
80       ! Compute residual, B_il
81       ! bottom interface l=1
82       !DIR$ SIMD
83       DO ij=ij_begin_ext,ij_end_ext
84          ml_g2 = gm2*m_il(ij,1) 
85          B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot
86          R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1))  &
87               + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-phis(ij)) )
88       ENDDO
89       ! inner interfaces
90       DO l=2,llm
91          !DIR$ SIMD
92          DO ij=ij_begin_ext,ij_end_ext
93             ml_g2 = gm2*m_il(ij,l) 
94             B_il(ij,l) = A_ik(ij,l)+A_ik(ij,l-1) + ml_g2
95             R_il(ij,l) = ml_g2*( Phi_il(ij,l)-Phi_star_il(ij,l))  &
96                     + tau2_g*(p_ik(ij,l)-p_ik(ij,l-1))
97             ! consistency check : if Wil=0 and initial state is in hydrostatic balance
98             ! then Phi_star_il(ij,l) = Phi_il(ij,l) - tau^2*g^2
99             ! and residual = tau^2*(ml+(1/g)dl_pi)=0
100          ENDDO
101       ENDDO
102       ! top interface l=llm+1
103       !DIR$ SIMD
104       DO ij=ij_begin_ext,ij_end_ext
105          ml_g2 = gm2*m_il(ij,llm+1) 
106          B_il(ij,llm+1) = A_ik(ij,llm) + ml_g2
107          R_il(ij,llm+1) = ml_g2*( Phi_il(ij,llm+1)-Phi_star_il(ij,llm+1))  &
108               + tau2_g*( ptop-p_ik(ij,llm) )
109       ENDDO
110
111       ! FIXME later
112       ! the lines below modify the tridiag problem
113       ! for flat, rigid boundary conditions at top and bottom :
114       ! zero out A(1), A(llm), R(1), R(llm+1)
115       ! => x(l)=0 at l=1,llm+1
116       DO ij=ij_begin_ext,ij_end_ext
117          A_ik(ij,1) = 0.
118          A_ik(ij,llm) = 0.
119          R_il(ij,1) = 0.
120          R_il(ij,llm+1) = 0.
121       ENDDO
122
123       IF(debug_hevi_solver) THEN ! print Linf(residual)
124          PRINT *, '[hevi_solver] R,p', iter, MAXVAL(ABS(R_il)), MAXVAL(p_ik)
125       END IF
126
127       ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
128       ! Forward sweep :
129       ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
130       ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
131       ! bottom interface l=1
132       !DIR$ SIMD
133       DO ij=ij_begin_ext,ij_end_ext
134          X_ij = 1./B_il(ij,1)
135          C_ik(ij,1) = -A_ik(ij,1) * X_ij 
136          D_il(ij,1) =  R_il(ij,1) * X_ij
137       ENDDO
138       ! inner interfaces/layers
139       DO l=2,llm
140          !DIR$ SIMD
141          DO ij=ij_begin_ext,ij_end_ext
142             X_ij = 1./(B_il(ij,l) + A_ik(ij,l-1)*C_ik(ij,l-1))
143             C_ik(ij,l) = -A_ik(ij,l) * X_ij 
144             D_il(ij,l) =  (R_il(ij,l)+A_ik(ij,l-1)*D_il(ij,l-1)) * X_ij
145          ENDDO
146       ENDDO
147       ! top interface l=llm+1
148       !DIR$ SIMD
149       DO ij=ij_begin_ext,ij_end_ext
150          X_ij = 1./(B_il(ij,llm+1) + A_ik(ij,llm)*C_ik(ij,llm))
151          D_il(ij,llm+1) =  (R_il(ij,llm+1)+A_ik(ij,llm)*D_il(ij,llm)) * X_ij
152       ENDDO
153       
154       ! Back substitution :
155       ! x(i) = D(i)-C(i)x(i+1), x(N+1)=0
156       ! + Newton-Raphson update
157       x_il=0. ! FIXME
158       ! top interface l=llm+1
159       !DIR$ SIMD
160       DO ij=ij_begin_ext,ij_end_ext
161          x_il(ij,llm+1) = D_il(ij,llm+1)
162          Phi_il(ij,llm+1) = Phi_il(ij,llm+1) - x_il(ij,llm+1)
163       ENDDO
164       ! lower interfaces
165       DO l=llm,1,-1
166          !DIR$ SIMD
167          DO ij=ij_begin_ext,ij_end_ext
168             x_il(ij,l) = D_il(ij,l) - C_ik(ij,l)*x_il(ij,l+1)
169             Phi_il(ij,l) = Phi_il(ij,l) - x_il(ij,l)
170          ENDDO
171       ENDDO
172
173       IF(debug_hevi_solver) THEN
174          PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
175          PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
176          DO l=1,llm+1
177             WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x', iter,l, MAXVAL(ABS(x_il(:,l)))
178          END DO
179       END IF
180
181    END DO ! Newton-Raphson
182
183    END IF ! dysl
184   
185  END SUBROUTINE compute_NH_geopot
186
187  SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)
188    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
189    REAL(rstd),INTENT(IN)    :: phis(iim*jjm)
190    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
191    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
192    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)
193    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
194    REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0
195    REAL(rstd),INTENT(OUT)   :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces
196    REAL(rstd),INTENT(OUT)   :: pres(iim*jjm,llm)          ! pressure
197    REAL(rstd),INTENT(OUT)   :: dW(iim*jjm,llm+1)
198    REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1)
199    REAL(rstd),INTENT(OUT)   :: du(3*iim*jjm,llm)
200
201    REAL(rstd) :: berni(iim*jjm,llm)         ! (W/m_il)^2
202    REAL(rstd) :: berni1(iim*jjm)         ! (W/m_il)^2
203    REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd, Rd_preff
204    INTEGER    :: ij, l
205
206    CALL trace_start("compute_caldyn_solver")
207
208    Rd=cpp*kappa
209
210    IF(dysl) THEN
211
212!$OMP BARRIER
213
214#include "../kernels_hex/caldyn_mil.k90"
215  IF(tau>0) THEN ! solve implicit problem for geopotential
216    CALL compute_NH_geopot(tau,phis, rhodz, m_il, theta, W, geopot)
217  END IF
218#define PHI_BOT(ij) phis(ij)
219#include "../kernels_hex/caldyn_solver.k90"
220#undef PHI_BOT
221!$OMP BARRIER
222
223    ELSE
224
225#define BERNI(ij) berni1(ij)
226    ! FIXME : vertical OpenMP parallelism will not work
227
228    ! average m_ik to interfaces => m_il
229    !DIR$ SIMD
230    DO ij=ij_begin_ext,ij_end_ext
231       m_il(ij,1) = .5*rhodz(ij,1)
232    ENDDO
233    DO l=2,llm
234       !DIR$ SIMD
235       DO ij=ij_begin_ext,ij_end_ext
236          m_il(ij,l) = .5*(rhodz(ij,l-1)+rhodz(ij,l))
237       ENDDO
238    ENDDO
239    !DIR$ SIMD
240    DO ij=ij_begin_ext,ij_end_ext
241       m_il(ij,llm+1) = .5*rhodz(ij,llm)
242    ENDDO
243
244    IF(tau>0) THEN ! solve implicit problem for geopotential
245       CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot)
246    END IF
247
248    ! Compute pressure, stored temporarily in pk
249    ! kappa = R/Cp
250    ! 1-kappa = Cv/Cp
251    ! Cp/Cv = 1/(1-kappa)
252    gamma = 1./(1.-kappa)
253    DO l=1,llm
254       !DIR$ SIMD
255       DO ij=ij_begin_ext,ij_end_ext
256          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))
257          X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij
258          ! kappa.theta.rho = p/exner
259          ! => X = (p/p0)/(exner/Cp)
260          !      = (p/p0)^(1-kappa)
261          pk(ij,l) = preff*(X_ij**gamma)
262       ENDDO
263    ENDDO
264
265    ! Update W, compute tendencies
266    DO l=2,llm
267       !DIR$ SIMD
268       DO ij=ij_begin_ext,ij_end_ext
269          dW(ij,l)   = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il(ij,l)
270          W(ij,l)    = W(ij,l)+tau*dW(ij,l) ! update W
271          dPhi(ij,l) = g*g*W(ij,l)/m_il(ij,l)
272       ENDDO
273       !       PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l)))
274       !       PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l)))
275    ENDDO
276    ! Lower BC (FIXME : no orography yet !)
277    DO ij=ij_begin,ij_end         
278       dPhi(ij,1)=0
279       W(ij,1)=0
280       dW(ij,1)=0
281       dPhi(ij,llm+1)=0 ! rigid lid
282       W(ij,llm+1)=0
283       dW(ij,llm+1)=0
284    ENDDO
285    ! Upper BC p=ptop
286    !    DO ij=ij_omp_begin_ext,ij_omp_end_ext
287    !       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm)
288    !       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm)
289    !    ENDDO
290   
291    ! Compute Exner function (needed by compute_caldyn_fast) and du=-g^2.grad(w^2)
292    DO l=1,llm
293       !DIR$ SIMD
294       DO ij=ij_begin_ext,ij_end_ext
295          pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow
296          BERNI(ij) = (-.25*g*g)*(              &
297                 (W(ij,l)/m_il(ij,l))**2       &
298               + (W(ij,l+1)/m_il(ij,l+1))**2 )
299       ENDDO
300       DO ij=ij_begin,ij_end         
301          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
302          du(ij+u_lup,l)   = ne_lup  *(BERNI(ij)-BERNI(ij+t_lup))
303          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
304       ENDDO
305    ENDDO
306#undef BERNI
307
308    END IF ! dysl
309
310    CALL trace_end("compute_caldyn_solver")
311   
312  END SUBROUTINE compute_caldyn_solver
313 
314  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du)
315    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
316    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0
317    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
318    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
319    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
320    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
321    REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm)
322    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
323    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
324
325    INTEGER :: i,j,ij,l
326    REAL(rstd) :: cp_ik, qv, temp, chi, nu, due, due_right, due_lup, due_ldown
327
328    CALL trace_start("compute_caldyn_fast")
329
330    IF(dysl_caldyn_fast) THEN
331#include "../kernels_hex/caldyn_fast.k90"
332    ELSE
333
334    ! Compute Bernoulli term
335    IF(boussinesq) THEN
336       DO l=ll_begin,ll_end
337          !DIR$ SIMD
338          DO ij=ij_begin,ij_end         
339             berni(ij,l) = pk(ij,l)
340             ! from now on pk contains the vertically-averaged geopotential
341             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
342          END DO
343       END DO
344    ELSE ! compressible
345
346       DO l=ll_begin,ll_end
347          SELECT CASE(caldyn_thermo)
348          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
349             !DIR$ SIMD
350             DO ij=ij_begin,ij_end         
351                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
352             END DO
353          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s)
354             !DIR$ SIMD
355             DO ij=ij_begin,ij_end         
356                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
357                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
358             END DO
359          CASE(thermo_moist) 
360             !DIR$ SIMD
361             DO ij=ij_begin,ij_end         
362                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
363                ! Bd = Phi + gibbs_d
364                ! Bv = Phi + gibbs_v
365                ! pk=temperature, theta=entropy
366                qv = theta(ij,l,2)
367                temp = pk(ij,l)
368                chi = log(temp/Treff)
369                nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
370                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
371                     + temp*(cpp*(1.-chi)+Rd*nu)
372                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
373                     + temp*(cppv*(1.-chi)+Rv*nu)
374            END DO
375          END SELECT
376       END DO
377
378    END IF ! Boussinesq/compressible
379
380!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
381    DO l=ll_begin,ll_end
382       IF(caldyn_thermo == thermo_moist) THEN
383          !DIR$ SIMD
384          DO ij=ij_begin,ij_end         
385             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
386                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
387                       *(pk(ij+t_right,l)-pk(ij,l))             &
388                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
389                       *(berniv(ij+t_right,l)-berniv(ij,l))
390             
391             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
392                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
393                       *(pk(ij+t_lup,l)-pk(ij,l))               &
394                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
395                       *(berniv(ij+t_lup,l)-berniv(ij,l))
396             
397             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
398                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
399                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
400                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
401                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
402             
403             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
404             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
405             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
406             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
407             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
408             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
409          END DO
410       ELSE
411          !DIR$ SIMD
412          DO ij=ij_begin,ij_end         
413             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
414                  *(pk(ij+t_right,l)-pk(ij,l))        &
415                  +  berni(ij+t_right,l)-berni(ij,l)
416             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
417                  *(pk(ij+t_lup,l)-pk(ij,l))          &
418                  +  berni(ij+t_lup,l)-berni(ij,l)
419             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
420                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
421                  +   berni(ij+t_ldown,l)-berni(ij,l)
422             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
423             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
424             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
425             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
426             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
427             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
428          END DO
429       END IF
430    END DO
431
432    END IF ! dysl
433    CALL trace_end("compute_caldyn_fast")
434
435  END SUBROUTINE compute_caldyn_fast
436
437  SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
438    REAL(rstd),INTENT(IN)  :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
439    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn) ! active scalars
440    REAL(rstd),INTENT(IN)  :: qu(3*iim*jjm,llm)
441    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
442    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
443    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
444   
445    REAL(rstd) :: Ftheta(3*iim*jjm,llm)  ! potential temperature flux
446    REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF
447    INTEGER :: ij,iq,l,kdown
448
449    CALL trace_start("compute_caldyn_Coriolis")
450
451    IF(dysl_caldyn_coriolis) THEN
452
453#include "../kernels_hex/coriolis.k90"
454
455    ELSE
456#define FTHETA(ij) Ftheta(ij,1)
457
458    DO l=ll_begin, ll_end
459       ! compute theta flux
460       DO iq=1,nqdyn
461       !DIR$ SIMD
462          DO ij=ij_begin_ext,ij_end_ext     
463             FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
464                                  * hflux(ij+u_right,l)
465             FTHETA(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
466                  * hflux(ij+u_lup,l)
467             FTHETA(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &
468                  * hflux(ij+u_ldown,l)
469          END DO
470          ! horizontal divergence of fluxes
471       !DIR$ SIMD
472          DO ij=ij_begin,ij_end         
473             ! dtheta_rhodz =  -div(flux.theta)
474             dtheta_rhodz(ij,l,iq)= &
475                  -1./Ai(ij)*(ne_right*FTHETA(ij+u_right)    +  &
476                  ne_rup*FTHETA(ij+u_rup)        +  & 
477                  ne_lup*FTHETA(ij+u_lup)        +  & 
478                  ne_left*FTHETA(ij+u_left)      +  & 
479                  ne_ldown*FTHETA(ij+u_ldown)    +  &
480                  ne_rdown*FTHETA(ij+u_rdown) )
481          END DO
482       END DO
483
484       !DIR$ SIMD
485       DO ij=ij_begin,ij_end         
486          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
487          convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
488               ne_rup*hflux(ij+u_rup,l)       +  & 
489               ne_lup*hflux(ij+u_lup,l)       +  & 
490               ne_left*hflux(ij+u_left,l)     +  & 
491               ne_ldown*hflux(ij+u_ldown,l)   +  & 
492               ne_rdown*hflux(ij+u_rdown,l))
493       END DO ! ij
494    END DO ! llm
495
496!!! Compute potential vorticity (Coriolis) contribution to du
497    SELECT CASE(caldyn_conserv)
498
499    CASE(conserv_energy) ! energy-conserving TRiSK
500
501       DO l=ll_begin,ll_end
502          !DIR$ SIMD
503          DO ij=ij_begin,ij_end         
504             uu_right = &
505                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
506                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
507                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
508                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
509                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
510                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_ldown,l))+         &
511                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rdown,l))+         &
512                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
513                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_rup,l))+             &
514                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_lup,l))
515             uu_lup = &
516                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
517                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
518                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
519                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
520                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
521                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_right,l)) +          &
522                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_rup,l)) +              &
523                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) +              &
524                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_left,l)) +            &
525                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_ldown,l))
526             uu_ldown = &
527                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
528                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
529                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
530                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
531                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
532                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_lup,l)) +          &
533                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_left,l)) +        &
534                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
535                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_rdown,l)) +      &
536                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_right,l))
537             du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right
538             du(ij+u_lup,l)   = du(ij+u_lup,l)   + .5*uu_lup
539             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + .5*uu_ldown
540          ENDDO
541       ENDDO
542
543    CASE(conserv_gassmann) ! energy-conserving TRiSK modified by Gassmann (2018)
544
545       DO l=ll_begin,ll_end
546          !DIR$ SIMD
547          DO ij=ij_begin,ij_end         
548             uu_right = &
549                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)  *qu(ij+t_right+u_lup,l)+                 &
550                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)  *qu(ij+u_rup,l)+                         &
551                .5*wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+     &
552                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*qu(ij+u_rdown,l)+                       &
553                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*qu(ij+t_right+u_ldown,l)+               &
554                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)*qu(ij+u_rdown,l)+               &
555                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)*qu(ij+t_right+u_ldown,l)+       &
556               .5*wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)*(qu(ij+u_right,l)+qu(ij+t_right+u_right,l))+         &
557                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)*qu(ij+t_right+u_lup,l)+           &
558                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)*qu(ij+u_rup,l)
559             uu_lup = &
560                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*qu(ij+t_lup+u_ldown,l) +                   &
561                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*qu(ij+u_left,l) +                         &
562               .5*wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +       &
563                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*qu(ij+u_rup,l) +                          &
564                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*qu(ij+t_lup+u_right,l)+                     & 
565                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)*qu(ij+u_rup,l) +                   &
566                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)*qu(ij+t_lup+u_right,l) +              &
567               .5*wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)*(qu(ij+u_lup,l)+qu(ij+t_lup+u_lup,l)) + &
568                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)*qu(ij+t_lup+u_ldown,l) +            &
569                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)*qu(ij+u_left,l)
570             uu_ldown = &
571                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*qu(ij+t_ldown,l+u_right) +               &
572                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*qu(ij+u_rdown,l) +                       &
573               .5*wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +        &
574                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*qu(ij+u_left,l) +                          &
575                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*qu(ij+t_ldown+u_lup,l) +                  & 
576                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)*qu(ij+u_left,l) +                    &
577                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)*qu(ij+t_ldown+u_lup,l) +          &
578               .5*wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)*(qu(ij+u_ldown,l)+qu(ij+t_ldown+u_ldown,l)) +      &
579                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)*qu(ij+t_ldown+u_right,l) +      &
580                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)*qu(ij+u_rdown,l)
581             du(ij+u_right,l) = du(ij+u_right,l) + uu_right
582             du(ij+u_lup,l)   = du(ij+u_lup,l)   + uu_lup
583             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + uu_ldown
584          ENDDO
585       ENDDO
586
587    CASE(conserv_enstrophy) ! enstrophy-conserving TRiSK
588
589       DO l=ll_begin,ll_end
590          !DIR$ SIMD
591          DO ij=ij_begin,ij_end         
592             uu_right = &
593                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
594                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
595                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
596                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
597                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
598                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
599                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
600                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
601                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
602                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
603             uu_lup = &
604                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
605                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
606                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
607                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
608                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
609                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
610                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
611                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
612                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
613                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
614             uu_ldown = &
615                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
616                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
617                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
618                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
619                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
620                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
621                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
622                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
623                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
624                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
625
626             du(ij+u_right,l) = du(ij+u_right,l) + uu_right*qu(ij+u_right,l)
627             du(ij+u_lup,l)   = du(ij+u_lup,l)   + uu_lup*qu(ij+u_lup,l)     
628             du(ij+u_ldown,l) = du(ij+u_ldown,l) + uu_ldown*qu(ij+u_ldown,l) 
629          END DO
630       END DO
631    CASE DEFAULT
632       STOP
633    END SELECT
634#undef FTHETA
635
636    END IF ! dysl
637
638    CALL trace_end("compute_caldyn_Coriolis")
639
640  END SUBROUTINE compute_caldyn_Coriolis
641
642  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hv, hflux,Kv,du, zero)
643    LOGICAL, INTENT(IN) :: zero
644    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
645    REAL(rstd),INTENT(IN)  :: Kv(2*iim*jjm,llm)   ! kinetic energy at vertices
646    REAL(rstd),INTENT(IN)  :: hv(2*iim*jjm,llm)   ! height/mass averaged to vertices
647    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
648    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
649    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
650   
651    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
652    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
653    REAL(rstd) :: uu_right, uu_lup, uu_ldown, ke, uu
654    INTEGER :: ij,l
655
656    CALL trace_start("compute_caldyn_slow_hydro")
657
658    IF(dysl_slow_hydro) THEN
659
660#define BERNI(ij,l) berni(ij,l)
661#include "../kernels_hex/caldyn_slow_hydro.k90"
662#undef BERNI
663
664     ELSE
665
666#define BERNI(ij) berni1(ij)
667
668    DO l = ll_begin, ll_end
669       !  Compute mass fluxes
670       IF (caldyn_conserv==conserv_energy) CALL test_message(req_qu) 
671
672       IF(caldyn_kinetic==kinetic_trisk) THEN
673          !DIR$ SIMD
674          DO ij=ij_begin_ext,ij_end_ext
675             uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)
676             uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)
677             uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)
678             uu_right= uu_right*le_de(ij+u_right)
679             uu_lup  = uu_lup  *le_de(ij+u_lup)
680             uu_ldown= uu_ldown*le_de(ij+u_ldown)
681             hflux(ij+u_right,l)=uu_right
682             hflux(ij+u_lup,l)  =uu_lup
683             hflux(ij+u_ldown,l)=uu_ldown
684          ENDDO
685       ELSE ! mass flux deriving from consistent kinetic energy
686          !DIR$ SIMD
687          DO ij=ij_begin_ext,ij_end_ext
688             uu_right=0.5*(hv(ij+z_rup,l)+hv(ij+z_rdown,l))*u(ij+u_right,l)
689             uu_lup=0.5*(hv(ij+z_up,l)+hv(ij+z_lup,l))*u(ij+u_lup,l)
690             uu_ldown=0.5*(hv(ij+z_ldown,l)+hv(ij+z_down,l))*u(ij+u_ldown,l)
691             uu_right= uu_right*le_de(ij+u_right)
692             uu_lup  = uu_lup  *le_de(ij+u_lup)
693             uu_ldown= uu_ldown*le_de(ij+u_ldown)
694             hflux(ij+u_right,l)=uu_right
695             hflux(ij+u_lup,l)  =uu_lup
696             hflux(ij+u_ldown,l)=uu_ldown
697          ENDDO
698       END IF
699
700       ! Compute Bernoulli=kinetic energy
701       IF(caldyn_kinetic==kinetic_trisk) THEN
702          !DIR$ SIMD
703          DO ij=ij_begin,ij_end         
704             BERNI(ij) = &
705                  1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
706                                le_de(ij+u_rup)*u(ij+u_rup,l)**2     +    &
707                                le_de(ij+u_lup)*u(ij+u_lup,l)**2     +    &
708                                le_de(ij+u_left)*u(ij+u_left,l)**2   +    &
709                                le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
710                                le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
711          ENDDO
712       ELSE
713          !DIR$ SIMD
714          DO ij=ij_begin,ij_end
715             BERNI(ij) = Riv(ij,vup)   *Kv(ij+z_up,l)    + &
716                         Riv(ij,vlup)  *Kv(ij+z_lup,l)   + &
717                         Riv(ij,vldown)*Kv(ij+z_ldown,l) + &
718                         Riv(ij,vdown) *Kv(ij+z_down,l)  + &
719                         Riv(ij,vrdown)*Kv(ij+z_rdown,l) + &
720                         Riv(ij,vrup)  *Kv(ij+z_rup,l)
721          END DO
722       END IF
723       ! Compute du=-grad(Bernoulli)
724       IF(zero) THEN
725          !DIR$ SIMD
726          DO ij=ij_begin,ij_end
727             du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
728             du(ij+u_lup,l)   = ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
729             du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
730          END DO
731       ELSE
732          !DIR$ SIMD
733          DO ij=ij_begin,ij_end
734             du(ij+u_right,l) = du(ij+u_right,l) + &
735                  ne_right*(BERNI(ij)-BERNI(ij+t_right))
736             du(ij+u_lup,l)   = du(ij+u_lup,l) + &
737                  ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
738             du(ij+u_ldown,l) = du(ij+u_ldown,l) + &
739                  ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
740          END DO
741       END IF
742    END DO
743
744#undef BERNI
745
746    END IF ! dysl
747    CALL trace_end("compute_caldyn_slow_hydro")   
748  END SUBROUTINE compute_caldyn_slow_hydro
749
750  SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
751    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
752    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)  ! rho*dz
753    REAL(rstd),INTENT(IN)  :: Phi(iim*jjm,llm+1)  ! prognostic geopotential
754    REAL(rstd),INTENT(IN)  :: W(iim*jjm,llm+1)    ! prognostic vertical momentum
755
756    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
757    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)
758    REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1)
759    REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1)
760
761    REAL(rstd) :: w_il(iim*jjm,llm+1) ! Wil/mil
762    REAL(rstd) :: F_el(3*iim*jjm,llm+1) ! NH mass flux
763    REAL(rstd) :: gradPhi2(iim*jjm,llm+1) ! grad_Phi**2
764    REAL(rstd) :: DePhil(3*iim*jjm,llm+1) ! grad(Phi)
765   
766    INTEGER :: ij,l,kdown,kup
767    REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu
768
769    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
770    REAL(rstd) :: G_el(3*iim*jjm,llm+1) ! horizontal flux of W
771    REAL(rstd) :: v_el(3*iim*jjm,llm+1)
772
773    REAL(rstd) :: berni1(iim*jjm)  ! Bernoulli function
774    REAL(rstd) :: G_el1(3*iim*jjm) ! horizontal flux of W
775    REAL(rstd) :: v_el1(3*iim*jjm)
776
777    CALL trace_start("compute_caldyn_slow_NH")
778
779    IF(dysl) THEN
780
781!$OMP BARRIER
782#include "../kernels_hex/caldyn_slow_NH.k90"
783!$OMP BARRIER
784     
785     ELSE
786
787#define BERNI(ij) berni1(ij)
788#define G_EL(ij) G_el1(ij)
789#define V_EL(ij) v_el1(ij)
790
791    DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
792       IF(l==1) THEN
793          kdown=1
794       ELSE
795          kdown=l-1
796       END IF
797       IF(l==llm+1) THEN
798          kup=llm
799       ELSE
800          kup=l
801       END IF
802       ! below : "checked" means "formula also valid when kup=kdown (top/bottom)"
803       ! compute mil, wil=Wil/mil
804       DO ij=ij_begin_ext, ij_end_ext
805          w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) ! checked
806       END DO
807       ! compute DePhi, v_el, G_el, F_el
808       ! v_el, W2_el and therefore G_el incorporate metric factor le_de
809       ! while DePhil, W_el and F_el don't
810       DO ij=ij_begin_ext, ij_end_ext
811          ! Compute on edge 'right'
812          W_el  = .5*( W(ij,l)+W(ij+t_right,l) )
813          DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
814          F_el(ij+u_right,l)   = DePhil(ij+u_right,l)*W_el
815          W2_el = .5*le_de(ij+u_right) * &
816               ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) )
817          V_EL(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked
818          G_EL(ij+u_right) = V_EL(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el
819          ! Compute on edge 'lup'
820          W_el  = .5*( W(ij,l)+W(ij+t_lup,l) )
821          DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
822          F_el(ij+u_lup,l)   = DePhil(ij+u_lup,l)*W_el
823          W2_el = .5*le_de(ij+u_lup) * &
824               ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) )
825          V_EL(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked
826          G_EL(ij+u_lup) = V_EL(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el
827          ! Compute on edge 'ldown'
828          W_el  = .5*( W(ij,l)+W(ij+t_ldown,l) )
829          DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
830          F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el
831          W2_el = .5*le_de(ij+u_ldown) * &
832               ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) )
833          V_EL(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked
834          G_EL(ij+u_ldown) = V_EL(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el
835       END DO
836       ! compute GradPhi2, dPhi, dW
837       DO ij=ij_begin_ext, ij_end_ext
838          gradPhi2(ij,l) = &
839               1/(2*Ai(ij))*(le_de(ij+u_right)*DePhil(ij+u_right,l)**2 + &
840               le_de(ij+u_rup)*DePhil(ij+u_rup,l)**2 +      &
841               le_de(ij+u_lup)*DePhil(ij+u_lup,l)**2 +      &
842               le_de(ij+u_left)*DePhil(ij+u_left,l)**2 +    &
843               le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 +  &
844               le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 )
845
846          dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* & 
847               ( DePhil(ij+u_right,l)*V_EL(ij+u_right) + & ! -v.gradPhi,
848                 DePhil(ij+u_rup,l)*V_EL(ij+u_rup) +     & ! v_el already has le_de
849                 DePhil(ij+u_lup,l)*V_EL(ij+u_lup) +     &
850                 DePhil(ij+u_left,l)*V_EL(ij+u_left) +   &
851                 DePhil(ij+u_ldown,l)*V_EL(ij+u_ldown) + &
852                 DePhil(ij+u_rdown,l)*V_EL(ij+u_rdown) )
853
854          dW(ij,l) = -1./Ai(ij)*(           & ! -div(G_el),
855               ne_right*G_EL(ij+u_right) +  & ! G_el already has le_de
856               ne_rup*G_EL(ij+u_rup) +      &
857               ne_lup*G_EL(ij+u_lup) +      & 
858               ne_left*G_EL(ij+u_left) +    &
859               ne_ldown*G_EL(ij+u_ldown) +  &
860               ne_rdown*G_EL(ij+u_rdown))
861       END DO
862    END DO
863
864    DO l=ll_begin, ll_end ! compute on k levels (layers)
865       ! Compute berni at scalar points
866       DO ij=ij_begin_ext, ij_end_ext
867          BERNI(ij) = &
868               1/(4*Ai(ij))*( &
869                    le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
870                    le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
871                    le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
872                    le_de(ij+u_left)*u(ij+u_left,l)**2 +       &
873                    le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
874                    le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) &
875               - .25*( gradPhi2(ij,l)  *w_il(ij,l)**2 +   &
876                      gradPhi2(ij,l+1)*w_il(ij,l+1)**2 )
877       END DO
878       ! Compute mass flux and grad(berni) at edges
879       DO ij=ij_begin_ext, ij_end_ext
880          ! Compute on edge 'right'
881          uu_right = 0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) &
882                    -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1))
883          hflux(ij+u_right,l) = uu_right*le_de(ij+u_right)
884          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
885          ! Compute on edge 'lup'
886          uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) &
887                  -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))
888          hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup)
889          du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup))
890          ! Compute on edge 'ldown'
891          uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) &
892                    -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1))
893          hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown)
894          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
895       END DO
896    END DO
897
898#undef V_EL
899#undef G_EL
900#undef BERNI
901
902    END IF ! dysl
903
904    CALL trace_end("compute_caldyn_slow_NH")
905
906  END SUBROUTINE compute_caldyn_slow_NH
907
908END MODULE caldyn_kernels_hevi_mod
Note: See TracBrowser for help on using the repository browser.