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

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

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