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

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

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