source: codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90 @ 734

Last change on this file since 734 was 734, checked in by dubos, 6 years ago

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