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

Last change on this file since 735 was 735, checked in by dubos, 6 years ago

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