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

Last change on this file since 844 was 844, checked in by jisesh, 5 years ago

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