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

Last change on this file since 404 was 404, checked in by dubos, 8 years ago

Towards moist dynamics - tested with DCMIP41
File size: 34.0 KB
Line 
1MODULE caldyn_kernels_hevi_mod
2  USE icosa
3  USE trace
4  USE omp_para
5  USE disvert_mod
6  USE transfert_mod
7  USE caldyn_kernels_base_mod
8  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) = (-.5*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)
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
390    INTEGER :: i,j,ij,l
391    REAL(rstd) :: due_right, due_lup, due_ldown
392
393    CALL trace_start("compute_caldyn_fast")
394
395    ! Compute Bernoulli term
396    IF(boussinesq) THEN
397
398       DO l=ll_begin,ll_end
399          !DIR$ SIMD
400          DO ij=ij_begin,ij_end         
401             berni(ij,l) = pk(ij,l)
402             ! from now on pk contains the vertically-averaged geopotential
403             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
404          END DO
405       END DO
406    ELSE ! compressible
407
408       DO l=ll_begin,ll_end
409          SELECT CASE(caldyn_thermo)
410          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
411             !DIR$ SIMD
412             DO ij=ij_begin,ij_end         
413                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
414             END DO
415          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s)
416             !DIR$ SIMD
417             DO ij=ij_begin,ij_end         
418                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
419                     + pk(ij,l)*(cpp-theta(ij,l)) ! pk=temperature, theta=entropy
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       !DIR$ SIMD
429       DO ij=ij_begin,ij_end         
430          due_right =  0.5*(theta(ij,l)+theta(ij+t_right,l))  &
431                          *(pk(ij+t_right,l)-pk(ij,l))        &
432                         +  berni(ij+t_right,l)-berni(ij,l)
433          due_lup = 0.5*(theta(ij,l)+theta(ij+t_lup,l))    &
434                       *(pk(ij+t_lup,l)-pk(ij,l))          &
435                         +  berni(ij+t_lup,l)-berni(ij,l)
436          due_ldown = 0.5*(theta(ij,l)+theta(ij+t_ldown,l)) &
437                         *(pk(ij+t_ldown,l)-pk(ij,l))      &
438                        +   berni(ij+t_ldown,l)-berni(ij,l)
439          du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
440          du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
441          du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
442          u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
443          u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
444          u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
445       ENDDO
446    ENDDO
447
448    CALL trace_end("compute_caldyn_fast")
449
450  END SUBROUTINE compute_caldyn_fast
451
452  SUBROUTINE compute_caldyn_Coriolis(hflux,theta,qu, convm,dtheta_rhodz,du)
453    REAL(rstd),INTENT(IN)  :: hflux(3*iim*jjm,llm)  ! hflux in kg/s
454    REAL(rstd),INTENT(IN)  :: theta(iim*jjm,llm,nqdyn) ! active scalars
455    REAL(rstd),INTENT(IN)  :: qu(3*iim*jjm,llm)
456    REAL(rstd),INTENT(OUT) :: convm(iim*jjm,llm)  ! mass flux convergence
457    REAL(rstd),INTENT(OUT) :: dtheta_rhodz(iim*jjm,llm,nqdyn)
458    REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm)
459   
460    REAL(rstd) :: Ftheta(3*iim*jjm)  ! potential temperature flux
461    REAL(rstd) :: uu_right, uu_lup, uu_ldown
462    INTEGER :: ij,iq,l,kdown
463
464    CALL trace_start("compute_caldyn_Coriolis")
465
466    DO l=ll_begin, ll_end
467       ! compute theta flux
468       !DIR$ SIMD
469       DO iq=1,nqdyn
470          DO ij=ij_begin_ext,ij_end_ext     
471             Ftheta(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &
472                                  * hflux(ij+u_right,l)
473             Ftheta(ij+u_lup)   = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &
474                  * hflux(ij+u_lup,l)
475             Ftheta(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &
476                  * hflux(ij+u_ldown,l)
477          END DO
478          ! horizontal divergence of fluxes
479          DO ij=ij_begin,ij_end         
480             ! dtheta_rhodz =  -div(flux.theta)
481             dtheta_rhodz(ij,l,iq)= &
482                  -1./Ai(ij)*(ne_right*Ftheta(ij+u_right)    +  &
483                  ne_rup*Ftheta(ij+u_rup)        +  & 
484                  ne_lup*Ftheta(ij+u_lup)        +  & 
485                  ne_left*Ftheta(ij+u_left)      +  & 
486                  ne_ldown*Ftheta(ij+u_ldown)    +  &
487                  ne_rdown*Ftheta(ij+u_rdown) )
488          END DO
489       END DO
490
491       DO ij=ij_begin,ij_end         
492          ! convm = -div(mass flux), sign convention as in Ringler et al. 2012, eq. 21
493          convm(ij,l)= -1./Ai(ij)*(ne_right*hflux(ij+u_right,l)   +  &
494               ne_rup*hflux(ij+u_rup,l)       +  & 
495               ne_lup*hflux(ij+u_lup,l)       +  & 
496               ne_left*hflux(ij+u_left,l)     +  & 
497               ne_ldown*hflux(ij+u_ldown,l)   +  & 
498               ne_rdown*hflux(ij+u_rdown,l))
499       END DO ! ij
500    END DO ! llm
501
502!!! Compute potential vorticity (Coriolis) contribution to du
503    SELECT CASE(caldyn_conserv)
504
505    CASE(energy) ! energy-conserving TRiSK
506
507       DO l=ll_begin,ll_end
508          !DIR$ SIMD
509          DO ij=ij_begin,ij_end         
510             uu_right = &
511                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)*(qu(ij+u_right,l)+qu(ij+u_rup,l))+                             &
512                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)*(qu(ij+u_right,l)+qu(ij+u_lup,l))+                             &
513                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)*(qu(ij+u_right,l)+qu(ij+u_left,l))+                           &
514                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)*(qu(ij+u_right,l)+qu(ij+u_ldown,l))+                         &
515                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)*(qu(ij+u_right,l)+qu(ij+u_rdown,l))+                         & 
516                  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))+         &
517                  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))+         &
518                  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))+         &
519                  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))+             &
520                  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))
521             uu_lup = &
522                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)*(qu(ij+u_lup,l)+qu(ij+u_left,l)) +                        &
523                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)*(qu(ij+u_lup,l)+qu(ij+u_ldown,l)) +                      &
524                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)*(qu(ij+u_lup,l)+qu(ij+u_rdown,l)) +                      &
525                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)*(qu(ij+u_lup,l)+qu(ij+u_right,l)) +                      &
526                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)*(qu(ij+u_lup,l)+qu(ij+u_rup,l)) +                          & 
527                  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)) +          &
528                  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)) +              &
529                  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)) +              &
530                  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)) +            &
531                  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))
532             uu_ldown = &
533                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)*(qu(ij+u_ldown,l)+qu(ij+u_rdown,l)) +                      &
534                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)*(qu(ij+u_ldown,l)+qu(ij+u_right,l)) +                      &
535                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)*(qu(ij+u_ldown,l)+qu(ij+u_rup,l)) +                          &
536                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)*(qu(ij+u_ldown,l)+qu(ij+u_lup,l)) +                          &
537                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)*(qu(ij+u_ldown,l)+qu(ij+u_left,l)) +                        & 
538                  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)) +          &
539                  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)) +        &
540                  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)) +      &
541                  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)) +      &
542                  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))
543             du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right
544             du(ij+u_lup,l)   = du(ij+u_lup,l)   + .5*uu_lup
545             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + .5*uu_ldown
546          ENDDO
547       ENDDO
548
549    CASE(enstrophy) ! enstrophy-conserving TRiSK
550
551       DO l=ll_begin,ll_end
552          !DIR$ SIMD
553          DO ij=ij_begin,ij_end         
554             uu_right = &
555                  wee(ij+u_right,1,1)*hflux(ij+u_rup,l)+                           &
556                  wee(ij+u_right,2,1)*hflux(ij+u_lup,l)+                           &
557                  wee(ij+u_right,3,1)*hflux(ij+u_left,l)+                          &
558                  wee(ij+u_right,4,1)*hflux(ij+u_ldown,l)+                         &
559                  wee(ij+u_right,5,1)*hflux(ij+u_rdown,l)+                         & 
560                  wee(ij+u_right,1,2)*hflux(ij+t_right+u_ldown,l)+                 &
561                  wee(ij+u_right,2,2)*hflux(ij+t_right+u_rdown,l)+                 &
562                  wee(ij+u_right,3,2)*hflux(ij+t_right+u_right,l)+                 &
563                  wee(ij+u_right,4,2)*hflux(ij+t_right+u_rup,l)+                   &
564                  wee(ij+u_right,5,2)*hflux(ij+t_right+u_lup,l)
565             uu_lup = &
566                  wee(ij+u_lup,1,1)*hflux(ij+u_left,l)+                        &
567                  wee(ij+u_lup,2,1)*hflux(ij+u_ldown,l)+                       &
568                  wee(ij+u_lup,3,1)*hflux(ij+u_rdown,l)+                       &
569                  wee(ij+u_lup,4,1)*hflux(ij+u_right,l)+                       &
570                  wee(ij+u_lup,5,1)*hflux(ij+u_rup,l)+                         & 
571                  wee(ij+u_lup,1,2)*hflux(ij+t_lup+u_right,l)+                 &
572                  wee(ij+u_lup,2,2)*hflux(ij+t_lup+u_rup,l)+                   &
573                  wee(ij+u_lup,3,2)*hflux(ij+t_lup+u_lup,l)+                   &
574                  wee(ij+u_lup,4,2)*hflux(ij+t_lup+u_left,l)+                  &
575                  wee(ij+u_lup,5,2)*hflux(ij+t_lup+u_ldown,l)
576             uu_ldown = &
577                  wee(ij+u_ldown,1,1)*hflux(ij+u_rdown,l)+                      &
578                  wee(ij+u_ldown,2,1)*hflux(ij+u_right,l)+                      &
579                  wee(ij+u_ldown,3,1)*hflux(ij+u_rup,l)+                        &
580                  wee(ij+u_ldown,4,1)*hflux(ij+u_lup,l)+                        &
581                  wee(ij+u_ldown,5,1)*hflux(ij+u_left,l)+                       & 
582                  wee(ij+u_ldown,1,2)*hflux(ij+t_ldown+u_lup,l)+                &
583                  wee(ij+u_ldown,2,2)*hflux(ij+t_ldown+u_left,l)+               &
584                  wee(ij+u_ldown,3,2)*hflux(ij+t_ldown+u_ldown,l)+              &
585                  wee(ij+u_ldown,4,2)*hflux(ij+t_ldown+u_rdown,l)+              &
586                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l)
587
588             du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right
589             du(ij+u_lup,l)   = du(ij+u_lup,l)   + .5*uu_lup
590             du(ij+u_ldown,l) = du(ij+u_ldown,l)   + .5*uu_ldown
591          END DO
592       END DO
593    CASE DEFAULT
594       STOP
595    END SELECT
596
597    CALL trace_end("compute_caldyn_Coriolis")
598
599  END SUBROUTINE compute_caldyn_Coriolis
600
601  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du)
602    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
603    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)
604    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
605    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)
606   
607    REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function
608    REAL(rstd) :: uu_right, uu_lup, uu_ldown
609    INTEGER :: ij,l
610
611    CALL trace_start("compute_caldyn_slow_hydro")
612
613    le_de(:) = le(:)/de(:) ! FIXME - make sure le_de is what we expect
614
615    DO l = ll_begin, ll_end
616       !  Compute mass fluxes
617       IF (caldyn_conserv==energy) CALL test_message(req_qu) 
618       !DIR$ SIMD
619       DO ij=ij_begin_ext,ij_end_ext
620          uu_right=0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l)
621          uu_lup=0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l)
622          uu_ldown=0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l)
623          uu_right= uu_right*le_de(ij+u_right)
624          uu_lup  = uu_lup  *le_de(ij+u_lup)
625          uu_ldown= uu_ldown*le_de(ij+u_ldown)
626          hflux(ij+u_right,l)=uu_right
627          hflux(ij+u_lup,l)  =uu_lup
628          hflux(ij+u_ldown,l)=uu_ldown
629       ENDDO
630       ! Compute Bernoulli=kinetic energy
631       !DIR$ SIMD
632       DO ij=ij_begin,ij_end         
633          berni(ij) = &
634               1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
635               le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
636               le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
637               le_de(ij+u_left)*u(ij+u_left,l)**2 +       &
638               le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
639               le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) 
640       ENDDO
641       ! Compute du=-grad(Bernoulli)
642       !DIR$ SIMD
643       DO ij=ij_begin,ij_end
644             du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right))
645             du(ij+u_lup,l)   = ne_lup*(berni(ij)-berni(ij+t_lup))
646             du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
647       END DO
648    END DO
649
650    CALL trace_end("compute_caldyn_slow_hydro")   
651  END SUBROUTINE compute_caldyn_slow_hydro
652
653  SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, hflux,du,dPhi,dW)
654    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity"
655    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm)  ! rho*dz
656    REAL(rstd),INTENT(IN)  :: Phi(iim*jjm,llm+1)  ! prognostic geopotential
657    REAL(rstd),INTENT(IN)  :: W(iim*jjm,llm+1)    ! prognostic vertical momentum
658
659    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s
660    REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)
661    REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1)
662    REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1)
663
664    REAL(rstd) :: w_il(3*iim*jjm,llm+1) ! Wil/mil
665    REAL(rstd) :: F_el(3*iim*jjm,llm+1) ! NH mass flux
666    REAL(rstd) :: GradPhi2(3*iim*jjm,llm+1) ! grad_Phi**2
667    REAL(rstd) :: DePhil(3*iim*jjm,llm+1) ! grad(Phi)
668    REAL(rstd) :: berni(iim*jjm)  ! Bernoulli function
669    REAL(rstd) :: G_el(3*iim*jjm) ! horizontal flux of W
670    REAL(rstd) :: v_el(3*iim*jjm)
671   
672    INTEGER :: ij,l,kdown,kup
673    REAL(rstd) :: uu_right, uu_lup, uu_ldown, W_el, W2_el
674
675    CALL trace_start("compute_caldyn_slow_NH")
676
677    le_de(:) = le(:)/de(:) ! FIXME - make sure le_de is what we expect
678
679    DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)
680       IF(l==1) THEN
681          kdown=1
682       ELSE
683          kdown=l-1
684       END IF
685       IF(l==llm+1) THEN
686          kup=llm
687       ELSE
688          kup=l
689       END IF
690       ! below : "checked" means "formula also valid when kup=kdown (top/bottom)"
691       ! compute mil, wil=Wil/mil
692       DO ij=ij_begin_ext, ij_end_ext
693          w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) ! checked
694       END DO
695       ! compute DePhi, v_el, G_el, F_el
696       ! v_el, W2_el and therefore G_el incorporate metric factor le_de
697       ! while DePhil, W_el and F_el don't
698       DO ij=ij_begin_ext, ij_end_ext
699          ! Compute on edge 'right'
700          W_el  = .5*( W(ij,l)+W(ij+t_right,l) )
701          DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))
702          F_el(ij+u_right,l)   = DePhil(ij+u_right,l)*W_el
703          W2_el = .5*le_de(ij+u_right) * &
704               ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) )
705          v_el(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked
706          G_el(ij+u_right) = v_el(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el
707          ! Compute on edge 'lup'
708          W_el  = .5*( W(ij,l)+W(ij+t_lup,l) )
709          DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))
710          F_el(ij+u_lup,l)   = DePhil(ij+u_lup,l)*W_el
711          W2_el = .5*le_de(ij+u_lup) * &
712               ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) )
713          v_el(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked
714          G_el(ij+u_lup) = v_el(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el
715          ! Compute on edge 'ldown'
716          W_el  = .5*( W(ij,l)+W(ij+t_ldown,l) )
717          DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))
718          F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el
719          W2_el = .5*le_de(ij+u_ldown) * &
720               ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) )
721          v_el(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked
722          G_el(ij+u_ldown) = v_el(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el
723       END DO
724       ! compute GradPhi2, dPhi, dW
725       DO ij=ij_begin_ext, ij_end_ext
726          gradPhi2(ij,l) = &
727               1/(2*Ai(ij))*(le_de(ij+u_right)*DePhil(ij+u_right,l)**2 + &
728               le_de(ij+u_rup)*DePhil(ij+u_rup,l)**2 +      &
729               le_de(ij+u_lup)*DePhil(ij+u_lup,l)**2 +      &
730               le_de(ij+u_left)*DePhil(ij+u_left,l)**2 +    &
731               le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 +  &
732               le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 )
733!          gradPhi2(ij,l) = 0. ! FIXME !!
734
735          dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* & 
736               ( DePhil(ij+u_right,l)*v_el(ij+u_right) + & ! -v.gradPhi,
737                 DePhil(ij+u_rup,l)*v_el(ij+u_rup) +     & ! v_el already has le_de
738                 DePhil(ij+u_lup,l)*v_el(ij+u_lup) +     &
739                 DePhil(ij+u_left,l)*v_el(ij+u_left) +   &
740                 DePhil(ij+u_ldown,l)*v_el(ij+u_ldown) + &
741                 DePhil(ij+u_rdown,l)*v_el(ij+u_rdown) )
742
743          dW(ij,l) = -1./Ai(ij)*(           & ! -div(G_el),
744               ne_right*G_el(ij+u_right) +  & ! G_el already has le_de
745               ne_rup*G_el(ij+u_rup) +      &
746               ne_lup*G_el(ij+u_lup) +      & 
747               ne_left*G_el(ij+u_left) +    &
748               ne_ldown*G_el(ij+u_ldown) +  &
749               ne_rdown*G_el(ij+u_rdown))
750       END DO
751    END DO
752    ! FIXME !!
753 !   F_el(:,:)=0.
754 !   dPhi(:,:)=0.
755 !   dW(:,:)=0.
756
757    DO l=ll_begin, ll_end ! compute on k levels (layers)
758       ! Compute berni at scalar points
759       DO ij=ij_begin_ext, ij_end_ext
760          berni(ij) = &
761               1/(4*Ai(ij))*( &
762                    le_de(ij+u_right)*u(ij+u_right,l)**2 +    &
763                    le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          &
764                    le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          &
765                    le_de(ij+u_left)*u(ij+u_left,l)**2 +       &
766                    le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    &
767                    le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) &
768               - .5*( gradPhi2(ij,l)  *w_il(ij,l)**2 +   &
769                      gradPhi2(ij,l+1)*w_il(ij,l+1)**2 )
770       END DO
771       ! Compute mass flux and grad(berni) at edges
772       DO ij=ij_begin_ext, ij_end_ext
773          ! Compute on edge 'right'
774          uu_right = 0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) &
775                    -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1))
776          hflux(ij+u_right,l) = uu_right*le_de(ij+u_right)
777          du(ij+u_right,l) = ne_right*(berni(ij)-berni(ij+t_right))
778          ! Compute on edge 'lup'
779          uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) &
780                  -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))
781          hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup)
782          du(ij+u_lup,l) = ne_lup*(berni(ij)-berni(ij+t_lup))
783          ! Compute on edge 'ldown'
784          uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) &
785                    -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1))
786          hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown)
787          du(ij+u_ldown,l) = ne_ldown*(berni(ij)-berni(ij+t_ldown))
788       END DO
789    END DO
790    ! FIXME !!
791    ! du(:,:)=0.
792    ! hflux(:,:)=0.
793
794    CALL trace_end("compute_caldyn_slow_NH")
795
796  END SUBROUTINE compute_caldyn_slow_NH
797
798END MODULE caldyn_kernels_hevi_mod
Note: See TracBrowser for help on using the repository browser.