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

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

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