source: codes/icosagcm/devel/src/kernels_unst/compute_NH_geopot.k90 @ 822

Last change on this file since 822 was 822, checked in by dubos, 5 years ago

devel : fix syntax of some WRITE statements (tolerated by ifort)
File size: 4.8 KB
Line 
1   !--------------------------------------------------------------------------
2   !---------------------------- compute_NH_geopot ----------------------------------
3   tau2_g=tau*tau/g
4   g2=g*g
5   gm2 = 1./g2
6   vreff = Treff*cpp/preff*kappa
7   gamma = 1./(1.-kappa)
8   !$OMP BARRIER
9   ! compute Phi_star
10   !$OMP DO SCHEDULE(STATIC)
11   DO ij=1,primal_num
12      DO l = 1,llm+1
13         Phi_star_il(l,ij) = Phi_il(l,ij) + tau*g2*(W_il(l,ij)/m_il(l,ij)-tau)
14      END DO
15   END DO
16   !$OMP END DO
17   ! Newton-Raphson iteration : Phi_il contains current guess value
18   DO iter=1,2 ! 2 iterations should be enough
19      ! Compute pressure, A_ik
20      SELECT CASE(caldyn_thermo)
21      CASE(thermo_theta)
22         !$OMP DO SCHEDULE(STATIC)
23         DO ij=1,primal_num
24            DO l = 1,llm
25               rho_ij = (g*m_ik(l,ij))/(Phi_il(l+1,ij)-Phi_il(l,ij))
26               X_ij = (cpp/preff)*kappa*theta(l,ij)*rho_ij
27               p_ik(l,ij) = preff*(X_ij**gamma)
28               c2_mik = gamma*p_ik(l,ij)/(rho_ij*m_ik(l,ij)) ! c^2 = gamma*R*T = gamma*p/rho
29               A_ik(l,ij) = c2_mik*(tau/g*rho_ij)**2
30            END DO
31         END DO
32         !$OMP END DO
33      CASE(thermo_entropy)
34         !$OMP DO SCHEDULE(STATIC)
35         DO ij=1,primal_num
36            DO l = 1,llm
37               rho_ij = (g*m_ik(l,ij))/(Phi_il(l+1,ij)-Phi_il(l,ij))
38               X_ij = log(vreff*rho_ij) + theta(l,ij)/cpp
39               p_ik(l,ij) = preff*exp(X_ij*gamma)
40               c2_mik = gamma*p_ik(l,ij)/(rho_ij*m_ik(l,ij)) ! c^2 = gamma*R*T = gamma*p/rho
41               A_ik(l,ij) = c2_mik*(tau/g*rho_ij)**2
42            END DO
43         END DO
44         !$OMP END DO
45      CASE DEFAULT
46         PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo
47         STOP
48      END SELECT
49      ! NB : A(1), A(llm), R(1), R(llm+1) = 0 => x(l)=0 at l=1,llm+1 => flat, rigid top and bottom
50      ! Solve -A(l-1)x(l-1) + B(l)x(l) - A(l)x(l+1) = R(l) using Thomas algorithm
51      !$OMP DO SCHEDULE(STATIC)
52      DO ij=1,primal_num
53         ! Compute residual R_il and B_il
54         ! bottom interface l=1
55         ml_g2 = gm2*m_il(1,ij)
56         B_il(1,ij) = A_ik(1,ij) + ml_g2 + tau2_g*rho_bot
57         R_il(1,ij) = ml_g2*( Phi_il(1,ij)-Phi_star_il(1,ij)) &
58         + tau2_g*( p_ik(1,ij)-pbot+rho_bot*(Phi_il(1,ij)-PHI_BOT(ij)) )
59         DO l = 2,llm
60            ! inner interfaces
61            ml_g2 = gm2*m_il(l,ij)
62            B_il(l,ij) = A_ik(l,ij)+A_ik(l-1,ij) + ml_g2
63            R_il(l,ij) = ml_g2*( Phi_il(l,ij)-Phi_star_il(l,ij)) &
64            + tau2_g*(p_ik(l,ij)-p_ik(l-1,ij))
65            ! consistency check : if Wil=0 and initial state is in hydrostatic balance
66            ! then Phi_star_il(l,ij) = Phi_il(l,ij) - tau^2*g^2
67            ! and residual = tau^2*(ml+(1/g)dl_pi)=0
68         END DO
69         ! top interface l=llm+1
70         ml_g2 = gm2*m_il(llm+1,ij)
71         B_il(llm+1,ij) = A_ik(llm+1 -1,ij) + ml_g2
72         R_il(llm+1,ij) = ml_g2*( Phi_il(llm+1,ij)-Phi_star_il(llm+1,ij)) &
73         + tau2_g*( ptop-p_ik(llm+1 -1,ij) )
74         !
75         ! Forward sweep :
76         ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
77         ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
78         X_ij = 1./B_il(1,ij)
79         C_ik(1,ij) = -A_ik(1,ij) * X_ij
80         D_il(1,ij) = R_il(1,ij) * X_ij
81         DO l = 2,llm
82            X_ij = 1./( B_il(l,ij) + A_ik(l-1,ij)*C_ik(l-1,ij) )
83            C_ik(l,ij) = -A_ik(l,ij) * X_ij
84            D_il(l,ij) = (R_il(l,ij)+A_ik(l-1,ij)*D_il(l-1,ij)) * X_ij
85         END DO
86         X_ij = 1./( B_il(llm+1,ij) + A_ik(llm+1 -1,ij)*C_ik(llm+1 -1,ij) )
87         D_il(llm+1,ij) = (R_il(llm+1,ij)+A_ik(llm+1 -1,ij)*D_il(llm+1 -1,ij)) * X_ij
88         ! Back substitution :
89         ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0
90         ! + Newton-Raphson update
91         ! top interface l=llm+1
92         x_il(llm+1,ij) = D_il(llm+1,ij)
93         Phi_il(llm+1,ij) = Phi_il(llm+1,ij) - x_il(llm+1,ij)
94         DO l = llm,1,-1
95            ! Back substitution at lower interfaces
96            x_il(l,ij) = D_il(l,ij) - C_ik(l,ij)*x_il(l+1,ij)
97            Phi_il(l,ij) = Phi_il(l,ij) - x_il(l,ij)
98         END DO
99      END DO
100      !$OMP END DO
101      IF(debug_hevi_solver) THEN
102         PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
103         PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
104         DO l=1,llm+1
105            WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x_il', iter,l, MAXVAL(ABS(x_il(l,:)))
106         END DO
107         DO l=1,llm+1
108            WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] R_il', iter,l, MAXVAL(ABS(R_il(l,:)))
109         END DO
110      END IF
111   END DO ! Newton-Raphson
112   !$OMP BARRIER
113   debug_hevi_solver=.FALSE.
114   !---------------------------- compute_NH_geopot ----------------------------------
115   !--------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.