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

Last change on this file since 536 was 536, checked in by dubos, 7 years ago

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