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

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

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