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

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

devel: separate module for compute_caldyn_solver and removed Rd duplication in it
File size: 11.5 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_fast,compute_NH_geopot
16
17CONTAINS
18
19  SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il)
20    REAL(rstd),INTENT(IN)    :: tau ! solve Phi-tau*dPhi/dt = Phi_rhs
21    REAL(rstd),INTENT(IN)    :: phis(iim*jjm)
22    REAL(rstd),INTENT(IN)    :: m_ik(iim*jjm,llm)
23    REAL(rstd),INTENT(IN)    :: m_il(iim*jjm,llm+1)
24    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm)
25    REAL(rstd),INTENT(IN)    :: W_il(iim*jjm,llm+1)   ! vertical momentum
26    REAL(rstd),INTENT(INOUT) :: Phi_il(iim*jjm,llm+1) ! geopotential
27
28    REAL(rstd) :: Phi_star_il(iim*jjm,llm+1) 
29    REAL(rstd) :: p_ik(iim*jjm,llm)      ! pressure
30    REAL(rstd) :: R_il(iim*jjm,llm+1)    ! rhs of tridiag problem
31    REAL(rstd) :: x_il(iim*jjm,llm+1)    ! solution of tridiag problem
32    REAL(rstd) :: A_ik(iim*jjm,llm)      ! off-diagonal coefficients of tridiag problem
33    REAL(rstd) :: B_il(iim*jjm,llm+1)    ! diagonal coefficients of tridiag problem
34    REAL(rstd) :: C_ik(iim*jjm,llm)      ! Thomas algorithm
35    REAL(rstd) :: D_il(iim*jjm,llm+1)    ! Thomas algorithm
36    REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij
37    REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik, vreff
38
39    INTEGER    :: iter, ij, l, ij_omp_begin_ext, ij_omp_end_ext
40
41    CALL distrib_level(ij_begin_ext,ij_end_ext, ij_omp_begin_ext,ij_omp_end_ext)
42
43    IF(dysl) THEN
44#define PHI_BOT(ij) phis(ij)
45#include "../kernels_hex/compute_NH_geopot.k90"
46#undef PHI_BOT
47    ELSE
48!    FIXME : vertical OpenMP parallelism will not work
49   
50    tau2_g=tau*tau/g
51    g2=g*g
52    gm2 = g**-2
53    gamma = 1./(1.-kappa)
54   
55    ! compute Phi_star
56    DO l=1,llm+1
57       !DIR$ SIMD
58       DO ij=ij_begin_ext,ij_end_ext
59          Phi_star_il(ij,l) = Phi_il(ij,l) + tau*g2*(W_il(ij,l)/m_il(ij,l)-tau)
60       ENDDO
61    ENDDO
62
63    ! Newton-Raphson iteration : Phi_il contains current guess value
64    DO iter=1,5 ! 2 iterations should be enough
65
66       ! Compute pressure, A_ik
67       DO l=1,llm
68          !DIR$ SIMD
69          DO ij=ij_begin_ext,ij_end_ext
70             rho_ij = (g*m_ik(ij,l))/(Phi_il(ij,l+1)-Phi_il(ij,l))
71             X_ij = (cpp/preff)*kappa*theta(ij,l)*rho_ij
72             p_ik(ij,l) = preff*(X_ij**gamma)
73             c2_mik = gamma*p_ik(ij,l)/(rho_ij*m_ik(ij,l)) ! c^2 = gamma*R*T = gamma*p/rho
74             A_ik(ij,l) = c2_mik*(tau/g*rho_ij)**2
75          ENDDO
76       ENDDO
77
78       ! Compute residual, B_il
79       ! bottom interface l=1
80       !DIR$ SIMD
81       DO ij=ij_begin_ext,ij_end_ext
82          ml_g2 = gm2*m_il(ij,1) 
83          B_il(ij,1) = A_ik(ij,1) + ml_g2 + tau2_g*rho_bot
84          R_il(ij,1) = ml_g2*( Phi_il(ij,1)-Phi_star_il(ij,1))  &
85               + tau2_g*( p_ik(ij,1)-pbot+rho_bot*(Phi_il(ij,1)-phis(ij)) )
86       ENDDO
87       ! inner interfaces
88       DO l=2,llm
89          !DIR$ SIMD
90          DO ij=ij_begin_ext,ij_end_ext
91             ml_g2 = gm2*m_il(ij,l) 
92             B_il(ij,l) = A_ik(ij,l)+A_ik(ij,l-1) + ml_g2
93             R_il(ij,l) = ml_g2*( Phi_il(ij,l)-Phi_star_il(ij,l))  &
94                     + tau2_g*(p_ik(ij,l)-p_ik(ij,l-1))
95             ! consistency check : if Wil=0 and initial state is in hydrostatic balance
96             ! then Phi_star_il(ij,l) = Phi_il(ij,l) - tau^2*g^2
97             ! and residual = tau^2*(ml+(1/g)dl_pi)=0
98          ENDDO
99       ENDDO
100       ! top interface l=llm+1
101       !DIR$ SIMD
102       DO ij=ij_begin_ext,ij_end_ext
103          ml_g2 = gm2*m_il(ij,llm+1) 
104          B_il(ij,llm+1) = A_ik(ij,llm) + ml_g2
105          R_il(ij,llm+1) = ml_g2*( Phi_il(ij,llm+1)-Phi_star_il(ij,llm+1))  &
106               + tau2_g*( ptop-p_ik(ij,llm) )
107       ENDDO
108
109       ! FIXME later
110       ! the lines below modify the tridiag problem
111       ! for flat, rigid boundary conditions at top and bottom :
112       ! zero out A(1), A(llm), R(1), R(llm+1)
113       ! => x(l)=0 at l=1,llm+1
114       DO ij=ij_begin_ext,ij_end_ext
115          A_ik(ij,1) = 0.
116          A_ik(ij,llm) = 0.
117          R_il(ij,1) = 0.
118          R_il(ij,llm+1) = 0.
119       ENDDO
120
121       IF(debug_hevi_solver) THEN ! print Linf(residual)
122          PRINT *, '[hevi_solver] R,p', iter, MAXVAL(ABS(R_il)), MAXVAL(p_ik)
123       END IF
124
125       ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
126       ! Forward sweep :
127       ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
128       ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
129       ! bottom interface l=1
130       !DIR$ SIMD
131       DO ij=ij_begin_ext,ij_end_ext
132          X_ij = 1./B_il(ij,1)
133          C_ik(ij,1) = -A_ik(ij,1) * X_ij 
134          D_il(ij,1) =  R_il(ij,1) * X_ij
135       ENDDO
136       ! inner interfaces/layers
137       DO l=2,llm
138          !DIR$ SIMD
139          DO ij=ij_begin_ext,ij_end_ext
140             X_ij = 1./(B_il(ij,l) + A_ik(ij,l-1)*C_ik(ij,l-1))
141             C_ik(ij,l) = -A_ik(ij,l) * X_ij 
142             D_il(ij,l) =  (R_il(ij,l)+A_ik(ij,l-1)*D_il(ij,l-1)) * X_ij
143          ENDDO
144       ENDDO
145       ! top interface l=llm+1
146       !DIR$ SIMD
147       DO ij=ij_begin_ext,ij_end_ext
148          X_ij = 1./(B_il(ij,llm+1) + A_ik(ij,llm)*C_ik(ij,llm))
149          D_il(ij,llm+1) =  (R_il(ij,llm+1)+A_ik(ij,llm)*D_il(ij,llm)) * X_ij
150       ENDDO
151       
152       ! Back substitution :
153       ! x(i) = D(i)-C(i)x(i+1), x(N+1)=0
154       ! + Newton-Raphson update
155       x_il=0. ! FIXME
156       ! top interface l=llm+1
157       !DIR$ SIMD
158       DO ij=ij_begin_ext,ij_end_ext
159          x_il(ij,llm+1) = D_il(ij,llm+1)
160          Phi_il(ij,llm+1) = Phi_il(ij,llm+1) - x_il(ij,llm+1)
161       ENDDO
162       ! lower interfaces
163       DO l=llm,1,-1
164          !DIR$ SIMD
165          DO ij=ij_begin_ext,ij_end_ext
166             x_il(ij,l) = D_il(ij,l) - C_ik(ij,l)*x_il(ij,l+1)
167             Phi_il(ij,l) = Phi_il(ij,l) - x_il(ij,l)
168          ENDDO
169       ENDDO
170
171       IF(debug_hevi_solver) THEN
172          PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
173          PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
174          DO l=1,llm+1
175             WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x', iter,l, MAXVAL(ABS(x_il(:,l)))
176          END DO
177       END IF
178
179    END DO ! Newton-Raphson
180
181    END IF ! dysl
182   
183  END SUBROUTINE compute_NH_geopot
184
185  SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du)
186    REAL(rstd),INTENT(IN)    :: tau                ! "solve" u-tau*du/dt = rhs
187    REAL(rstd),INTENT(INOUT) :: u(iim*3*jjm,llm)   ! OUT if tau>0
188    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
189    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
190    REAL(rstd),INTENT(INOUT) :: pk(iim*jjm,llm)
191    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
192    REAL(rstd),INTENT(INOUT)   :: du(iim*3*jjm,llm)
193    REAL(rstd) :: berni(iim*jjm,llm)  ! Bernoulli function
194    REAL(rstd) :: berniv(iim*jjm,llm)  ! moist Bernoulli function
195
196    INTEGER :: i,j,ij,l
197    REAL(rstd) :: cp_ik, qv, temp, chi, nu, due, due_right, due_lup, due_ldown
198
199    CALL trace_start("compute_caldyn_fast")
200
201    IF(dysl_caldyn_fast) THEN
202#include "../kernels_hex/caldyn_fast.k90"
203    ELSE
204
205    ! Compute Bernoulli term
206    IF(boussinesq) THEN
207       DO l=ll_begin,ll_end
208          !DIR$ SIMD
209          DO ij=ij_begin,ij_end         
210             berni(ij,l) = pk(ij,l)
211             ! from now on pk contains the vertically-averaged geopotential
212             pk(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
213          END DO
214       END DO
215    ELSE ! compressible
216
217       DO l=ll_begin,ll_end
218          SELECT CASE(caldyn_thermo)
219          CASE(thermo_theta) ! vdp = theta.dpi => B = Phi
220             !DIR$ SIMD
221             DO ij=ij_begin,ij_end         
222                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1))
223             END DO
224          CASE(thermo_entropy) ! vdp = dG + sdT => B = Phi + G, G=h-Ts=T*(cpp-s)
225             !DIR$ SIMD
226             DO ij=ij_begin,ij_end         
227                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
228                     + pk(ij,l)*(cpp-theta(ij,l,1)) ! pk=temperature, theta=entropy
229             END DO
230          CASE(thermo_moist) 
231             !DIR$ SIMD
232             DO ij=ij_begin,ij_end         
233                ! du/dt = grad(Bd)+rv.grad(Bv)+s.grad(T)
234                ! Bd = Phi + gibbs_d
235                ! Bv = Phi + gibbs_v
236                ! pk=temperature, theta=entropy
237                qv = theta(ij,l,2)
238                temp = pk(ij,l)
239                chi = log(temp/Treff)
240                nu = (chi*(cpp+qv*cppv)-theta(ij,l,1))/(Rd+qv*Rv) ! log(p/preff)
241                berni(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
242                     + temp*(cpp*(1.-chi)+Rd*nu)
243                berniv(ij,l) = .5*(geopot(ij,l)+geopot(ij,l+1)) &
244                     + temp*(cppv*(1.-chi)+Rv*nu)
245            END DO
246          END SELECT
247       END DO
248
249    END IF ! Boussinesq/compressible
250
251!!! u:=u+tau*du, du = -grad(B)-theta.grad(pi)
252    DO l=ll_begin,ll_end
253       IF(caldyn_thermo == thermo_moist) THEN
254          !DIR$ SIMD
255          DO ij=ij_begin,ij_end         
256             due_right =  berni(ij+t_right,l)-berni(ij,l)       &
257                  + 0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))   &
258                       *(pk(ij+t_right,l)-pk(ij,l))             &
259                  + 0.5*(theta(ij,l,2)+theta(ij+t_right,l,2))   &
260                       *(berniv(ij+t_right,l)-berniv(ij,l))
261             
262             due_lup = berni(ij+t_lup,l)-berni(ij,l)            &
263                  + 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))     &
264                       *(pk(ij+t_lup,l)-pk(ij,l))               &
265                  + 0.5*(theta(ij,l,2)+theta(ij+t_lup,l,2))     &
266                       *(berniv(ij+t_lup,l)-berniv(ij,l))
267             
268             due_ldown = berni(ij+t_ldown,l)-berni(ij,l)        &
269                  + 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1))   &
270                       *(pk(ij+t_ldown,l)-pk(ij,l))             &
271                  + 0.5*(theta(ij,l,2)+theta(ij+t_ldown,l,2))   &
272                       *(berniv(ij+t_ldown,l)-berniv(ij,l))
273             
274             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
275             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
276             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
277             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
278             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
279             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
280          END DO
281       ELSE
282          !DIR$ SIMD
283          DO ij=ij_begin,ij_end         
284             due_right =  0.5*(theta(ij,l,1)+theta(ij+t_right,l,1))  &
285                  *(pk(ij+t_right,l)-pk(ij,l))        &
286                  +  berni(ij+t_right,l)-berni(ij,l)
287             due_lup = 0.5*(theta(ij,l,1)+theta(ij+t_lup,l,1))    &
288                  *(pk(ij+t_lup,l)-pk(ij,l))          &
289                  +  berni(ij+t_lup,l)-berni(ij,l)
290             due_ldown = 0.5*(theta(ij,l,1)+theta(ij+t_ldown,l,1)) &
291                  *(pk(ij+t_ldown,l)-pk(ij,l))      &
292                  +   berni(ij+t_ldown,l)-berni(ij,l)
293             du(ij+u_right,l) = du(ij+u_right,l) - ne_right*due_right
294             du(ij+u_lup,l)   = du(ij+u_lup,l)   - ne_lup*due_lup
295             du(ij+u_ldown,l) = du(ij+u_ldown,l) - ne_ldown*due_ldown
296             u(ij+u_right,l) = u(ij+u_right,l) + tau*du(ij+u_right,l)
297             u(ij+u_lup,l)   = u(ij+u_lup,l)   + tau*du(ij+u_lup,l)
298             u(ij+u_ldown,l) = u(ij+u_ldown,l) + tau*du(ij+u_ldown,l)
299          END DO
300       END IF
301    END DO
302
303    END IF ! dysl
304    CALL trace_end("compute_caldyn_fast")
305
306  END SUBROUTINE compute_caldyn_fast
307
308END MODULE caldyn_kernels_hevi_mod
Note: See TracBrowser for help on using the repository browser.