source: codes/icosagcm/devel/src/dynamics/compute_NH_geopot.F90 @ 915

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

devel: moved DYSL into compute_caldyn_solver.F90,compute_theta.F90 and compute_NH_geopot.F90
File size: 11.8 KB
Line 
1MODULE compute_NH_geopot_mod
2  USE grid_param, ONLY : llm
3  IMPLICIT NONE
4  PRIVATE
5
6  LOGICAL, SAVE :: debug_hevi_solver = .FALSE.
7
8#include "../unstructured/unstructured.h90"
9
10  PUBLIC :: compute_NH_geopot,compute_NH_geopot_unst
11
12CONTAINS
13
14#ifdef BEGIN_DYSL
15
16KERNEL(compute_NH_geopot)
17  tau2_g=tau*tau/g
18  g2=g*g
19  gm2 = 1./g2
20  vreff = Treff*cpp/preff*kappa
21  gamma = 1./(1.-kappa)
22
23  BARRIER
24
25  ! compute Phi_star
26  SEQUENCE_C1
27    BODY('1,llm+1')
28      Phi_star_il(CELL) = Phi_il(CELL) + tau*g2*(W_il(CELL)/m_il(CELL)-tau)
29    END_BLOCK
30  END_BLOCK
31
32  ! Newton-Raphson iteration : Phi_il contains current guess value
33  DO iter=1,2 ! 2 iterations should be enough
34    ! Compute pressure, A_ik
35    SELECT CASE(caldyn_thermo)
36    CASE(thermo_theta)
37      {% call() compute_p_and_Aik() %}
38        X_ij = (cpp/preff)*kappa*theta(CELL)*rho_ij
39        p_ik(CELL) = preff*(X_ij**gamma)
40        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho
41      {% endcall %}
42    CASE(thermo_entropy)
43      {% call() compute_p_and_Aik() %}
44        X_ij = log(vreff*rho_ij) + theta(CELL)/cpp
45        p_ik(CELL) = preff*exp(X_ij*gamma)
46        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho
47      {% endcall %}
48    CASE DEFAULT
49        PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo
50        STOP
51    END SELECT
52
53    ! NB : A(1), A(llm), R(1), R(llm+1) = 0 => x(l)=0 at l=1,llm+1 => flat, rigid top and bottom
54    ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
55
56    SEQUENCE_C1
57      ! Compute residual R_il and B_il
58      PROLOGUE(1)
59        ! bottom interface l=1
60        ml_g2 = gm2*m_il(CELL)
61        B_il(CELL) = A_ik(CELL) + ml_g2 + tau2_g*rho_bot
62        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
63                  + tau2_g*( p_ik(CELL)-pbot+rho_bot*(Phi_il(CELL)-PHI_BOT(HIDX(CELL))) )
64      END_BLOCK
65      BODY('2,llm')
66        ! inner interfaces
67        ml_g2 = gm2*m_il(CELL)
68        B_il(CELL) = A_ik(CELL)+A_ik(DOWN(CELL)) + ml_g2
69        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
70                  + tau2_g*(p_ik(CELL)-p_ik(DOWN(CELL)))
71        ! consistency check : if Wil=0 and initial state is in hydrostatic balance
72        ! then Phi_star_il(CELL) = Phi_il(CELL) - tau^2*g^2
73        ! and residual = tau^2*(ml+(1/g)dl_pi)=0
74      END_BLOCK
75      EPILOGUE('llm+1')
76        ! top interface l=llm+1
77        ml_g2 = gm2*m_il(CELL)
78        B_il(CELL) = A_ik(DOWN(CELL)) + ml_g2
79        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
80                  + tau2_g*( ptop-p_ik(DOWN(CELL)) )
81      END_BLOCK
82      !
83      ! Forward sweep :
84      ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
85      ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
86      PROLOGUE(1)
87        X_ij = 1./B_il(CELL)
88        C_ik(CELL) = -A_ik(CELL) * X_ij
89        D_il(CELL) =  R_il(CELL) * X_ij
90      END_BLOCK
91      BODY('2,llm')
92        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )
93        C_ik(CELL) = -A_ik(CELL) * X_ij
94        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij
95      END_BLOCK
96      EPILOGUE('llm+1')
97        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )
98        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij
99        ! Back substitution :
100        ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0
101        ! + Newton-Raphson update
102        ! top interface l=llm+1
103        x_il(CELL) = D_il(CELL)
104        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)
105      END_BLOCK
106      BODY('llm,1,-1')
107        ! Back substitution at lower interfaces
108        x_il(CELL) = D_il(CELL) - C_ik(CELL)*x_il(UP(CELL))
109        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)
110      END_BLOCK
111    END_BLOCK
112
113    IF(debug_hevi_solver) THEN
114       PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
115       PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
116       DO l=1,llm+1
117          WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x_il', iter,l, MAXVAL(ABS(x_il(l,:)))
118       END DO
119       DO l=1,llm+1
120          WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] R_il', iter,l, MAXVAL(ABS(R_il(l,:)))
121       END DO
122    END IF
123
124  END DO ! Newton-Raphson
125
126  BARRIER
127
128  debug_hevi_solver=.FALSE.
129END_BLOCK
130
131#endif END_DYSL
132
133SUBROUTINE compute_NH_geopot_unst(tau, m_ik, m_il, theta, W_il, Phi_il)
134  USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
135  USE disvert_mod, only : g,Treff,cpp,preff,kappa,caldyn_thermo,thermo_theta, &
136    thermo_entropy
137  USE disvert_mod, ONLY : ptop
138  USE data_unstructured_mod, ONLY : primal_num,edge_num,dual_num,id_NH_geopot,debug_hevi_solver_, &
139    PHI_BOT,pbot,rho_bot,enter_trace, exit_trace
140  FIELD_MASS   :: m_ik, theta   ! IN*2
141  FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il  ! IN,INOUT*2, LOCAL
142  NUM :: tau, gamma, tau2_g, tau2_g2, g2, gm2, vreff, Rd_preff
143  INTEGER :: iter
144  LOGICAL :: debug_hevi_solver
145  DECLARE_INDICES
146  NUM :: rho_ij, X_ij, Y_ij, wil, rho_c2_mik, c2_mik, ml_g2
147#define COLUMN 0
148#if COLUMN 
149  NUM1(llm)  :: pk, Ak, Ck
150  NUM1(llm+1):: Rl, Bl, Dl, xl
151#define p_ik(l,ij) pk(l)
152#define A_ik(l,ij) Ak(l)
153#define C_ik(l,ij) Ck(l)
154#define R_il(l,ij) Rl(l)
155#define B_il(l,ij) Bl(l)
156#define D_il(l,ij) Dl(l)
157#define x_il(l,ij) xl(l)
158#else
159  FIELD_MASS :: p_ik, A_ik, C_ik
160  FIELD_GEOPOT :: R_il, B_il, D_il, x_il
161#endif
162
163  debug_hevi_solver=.FALSE.
164!$OMP MASTER
165  debug_hevi_solver = debug_hevi_solver_
166!$OMP END MASTER
167
168#define PHI_BOT(ij) Phi_bot
169  START_TRACE(id_NH_geopot, 7,0,0)
170#include "../kernels_unst/compute_NH_geopot.k90"
171  STOP_TRACE
172
173!$OMP MASTER
174  debug_hevi_solver_ = debug_hevi_solver
175!$OMP END MASTER
176#undef PHI_BOT
177
178#if COLUMN
179#undef p_ik
180#undef A_ik
181#undef C_ik
182#undef R_il
183#undef B_il
184#undef D_il
185#undef x_il
186#endif
187#undef COLUMN
188END SUBROUTINE compute_NH_geopot_unst
189
190  SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il)
191    USE icosa
192    USE disvert_mod
193    USE caldyn_vars_mod
194    USE omp_para
195    REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6
196    REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs
197    REAL(rstd),INTENT(IN)    :: phis(iim*jjm)
198    REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm)
199    REAL(rstd),INTENT(IN)    :: m_il(iim*jjm,llm+1)
200    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
201    REAL(rstd),INTENT(IN)    :: W_il(iim*jjm,llm+1)   ! vertical momentum
202    REAL(rstd),INTENT(INOUT) :: Phi_il(iim*jjm,llm+1) ! geopotential
203
204    REAL(rstd) :: Phi_star_il(iim*jjm,llm+1) 
205    REAL(rstd) :: p_ik(iim*jjm,llm)      ! pressure
206    REAL(rstd) :: R_il(iim*jjm,llm+1)    ! rhs of tridiag problem
207    REAL(rstd) :: x_il(iim*jjm,llm+1)    ! solution of tridiag problem
208    REAL(rstd) :: A_ik(iim*jjm,llm)      ! off-diagonal coefficients of tridiag problem
209    REAL(rstd) :: B_il(iim*jjm,llm+1)    ! diagonal coefficients of tridiag problem
210    REAL(rstd) :: C_ik(iim*jjm,llm)      ! Thomas algorithm
211    REAL(rstd) :: D_il(iim*jjm,llm+1)    ! Thomas algorithm
212    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij
213    REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik, vreff
214
215    INTEGER    :: iter, ij, l, ij_omp_begin_ext, ij_omp_end_ext
216
217    CALL distrib_level(ij_begin_ext,ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext)
218
219    IF(dysl) THEN
220#define PHI_BOT(ij) phis(ij)
221#include "../kernels_hex/compute_NH_geopot.k90"
222#undef PHI_BOT
223    ELSE
224!    FIXME : vertical OpenMP parallelism will not work
225   
226    tau2_g=tau*tau/g
227    g2=g*g
228    gm2 = g**-2
229    gamma = 1./(1.-kappa)
230   
231    ! compute Phi_star
232    DO l=1,llm+1
233       !DIR$ SIMD
234       DO ij=ij_begin_ext,ij_end_ext
235          Phi_star_il(ij,l) = Phi_il(ij,l) + tau*g2*(W_il(ij,l)/m_il(ij,l)-tau)
236       ENDDO
237    ENDDO
238
239    ! Newton-Raphson iteration : Phi_il contains current guess value
240    DO iter=1,5 ! 2 iterations should be enough
241
242       ! Compute pressure, A_ik
243       DO l=1,llm
244          !DIR$ SIMD
245          DO ij=ij_begin_ext,ij_end_ext
246             rho_ij = (g*m_ik(ij,l))/(Phi_il(ij,l+1)-Phi_il(ij,l))
247             X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij
248             p_ik(ij,l) = preff*(X_ij**gamma)
249             c2_mik = gamma*p_ik(ij,l)/(rho_ij*m_ik(ij,l)) ! c^2 = gamma*R*T = gamma*p/rho
250             A_ik(ij,l) = c2_mik*(tau/g*rho_ij)**2
251          ENDDO
252       ENDDO
253
254       ! Compute residual, B_il
255       ! bottom interface l=1
256       !DIR$ SIMD
257       DO ij=ij_begin_ext,ij_end_ext
258          ml_g2 = gm2*m_il(ij,1) 
259          B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot
260          R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1))  &
261               + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-phis(ij)) )
262       ENDDO
263       ! inner interfaces
264       DO l=2,llm
265          !DIR$ SIMD
266          DO ij=ij_begin_ext,ij_end_ext
267             ml_g2 = gm2*m_il(ij,l) 
268             B_il(ij,l) = A_ik(ij,l)+A_ik(ij,l-1) + ml_g2
269             R_il(ij,l) = ml_g2*( Phi_il(ij,l)-Phi_star_il(ij,l))  &
270                     + tau2_g*(p_ik(ij,l)-p_ik(ij,l-1))
271             ! consistency check : if Wil=0 and initial state is in hydrostatic balance
272             ! then Phi_star_il(ij,l) = Phi_il(ij,l) - tau^2*g^2
273             ! and residual = tau^2*(ml+(1/g)dl_pi)=0
274          ENDDO
275       ENDDO
276       ! top interface l=llm+1
277       !DIR$ SIMD
278       DO ij=ij_begin_ext,ij_end_ext
279          ml_g2 = gm2*m_il(ij,llm+1) 
280          B_il(ij,llm+1) = A_ik(ij,llm) + ml_g2
281          R_il(ij,llm+1) = ml_g2*( Phi_il(ij,llm+1)-Phi_star_il(ij,llm+1))  &
282               + tau2_g*( ptop-p_ik(ij,llm) )
283       ENDDO
284
285       ! FIXME later
286       ! the lines below modify the tridiag problem
287       ! for flat, rigid boundary conditions at top and bottom :
288       ! zero out A(1), A(llm), R(1), R(llm+1)
289       ! => x(l)=0 at l=1,llm+1
290       DO ij=ij_begin_ext,ij_end_ext
291          A_ik(ij,1) = 0.
292          A_ik(ij,llm) = 0.
293          R_il(ij,1) = 0.
294          R_il(ij,llm+1) = 0.
295       ENDDO
296
297       IF(debug_hevi_solver) THEN ! print Linf(residual)
298          PRINT *, '[hevi_solver] R,p', iter, MAXVAL(ABS(R_il)), MAXVAL(p_ik)
299       END IF
300
301       ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
302       ! Forward sweep :
303       ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
304       ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
305       ! bottom interface l=1
306       !DIR$ SIMD
307       DO ij=ij_begin_ext,ij_end_ext
308          X_ij = 1./B_il(ij,1)
309          C_ik(ij,1) = -A_ik(ij,1) * X_ij 
310          D_il(ij,1) =  R_il(ij,1) * X_ij
311       ENDDO
312       ! inner interfaces/layers
313       DO l=2,llm
314          !DIR$ SIMD
315          DO ij=ij_begin_ext,ij_end_ext
316             X_ij = 1./(B_il(ij,l) + A_ik(ij,l-1)*C_ik(ij,l-1))
317             C_ik(ij,l) = -A_ik(ij,l) * X_ij 
318             D_il(ij,l) =  (R_il(ij,l)+A_ik(ij,l-1)*D_il(ij,l-1)) * X_ij
319          ENDDO
320       ENDDO
321       ! top interface l=llm+1
322       !DIR$ SIMD
323       DO ij=ij_begin_ext,ij_end_ext
324          X_ij = 1./(B_il(ij,llm+1) + A_ik(ij,llm)*C_ik(ij,llm))
325          D_il(ij,llm+1) =  (R_il(ij,llm+1)+A_ik(ij,llm)*D_il(ij,llm)) * X_ij
326       ENDDO
327       
328       ! Back substitution :
329       ! x(i) = D(i)-C(i)x(i+1), x(N+1)=0
330       ! + Newton-Raphson update
331       x_il=0. ! FIXME
332       ! top interface l=llm+1
333       !DIR$ SIMD
334       DO ij=ij_begin_ext,ij_end_ext
335          x_il(ij,llm+1) = D_il(ij,llm+1)
336          Phi_il(ij,llm+1) = Phi_il(ij,llm+1) - x_il(ij,llm+1)
337       ENDDO
338       ! lower interfaces
339       DO l=llm,1,-1
340          !DIR$ SIMD
341          DO ij=ij_begin_ext,ij_end_ext
342             x_il(ij,l) = D_il(ij,l) - C_ik(ij,l)*x_il(ij,l+1)
343             Phi_il(ij,l) = Phi_il(ij,l) - x_il(ij,l)
344          ENDDO
345       ENDDO
346
347       IF(debug_hevi_solver) THEN
348          PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
349          PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
350          DO l=1,llm+1
351             WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x', iter,l, MAXVAL(ABS(x_il(:,l)))
352          END DO
353       END IF
354
355    END DO ! Newton-Raphson
356
357    END IF ! dysl
358   
359  END SUBROUTINE compute_NH_geopot
360
361END MODULE compute_NH_geopot_mod
Note: See TracBrowser for help on using the repository browser.