source: codes/icosagcm/trunk/src/caldyn_kernels_hevi.f90 @ 519

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

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