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

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

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