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

Last change on this file since 1026 was 1026, checked in by dubos, 4 years ago

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