source: codes/icosagcm/trunk/src/dynamics/caldyn_kernels_hevi.F90

Last change on this file was 954, checked in by adurocher, 5 years ago

trunk : Added metric terms to kernels parameters to avoid Host/GPU transferts

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