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

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

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