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

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

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