source: codes/icosagcm/devel/Python/src/kernels_caldyn_NH.jin @ 876

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

devel: moved DYSL into compute_caldyn_slow_NH.F90 and compute_caldyn_Coriolis.F90
File size: 8.2 KB
Line 
1{% macro compute_p_and_Aik() %}
2{% set p_and_c2_from_theta=caller %}
3    SEQUENCE_C1
4      BODY('1,llm')
5          rho_ij = (g*m_ik(CELL))/(Phi_il(UP(CELL))-Phi_il(CELL))
6          {{ p_and_c2_from_theta() }}
7          A_ik(CELL) = c2_mik*(tau/g*rho_ij)**2
8      END_BLOCK
9    END_BLOCK
10{%- endmacro %}
11
12KERNEL(compute_NH_geopot)
13  tau2_g=tau*tau/g
14  g2=g*g
15  gm2 = 1./g2
16  vreff = Treff*cpp/preff*kappa
17  gamma = 1./(1.-kappa)
18   
19  BARRIER
20
21  ! compute Phi_star
22  SEQUENCE_C1
23    BODY('1,llm+1')
24      Phi_star_il(CELL) = Phi_il(CELL) + tau*g2*(W_il(CELL)/m_il(CELL)-tau)
25    END_BLOCK
26  END_BLOCK
27
28  ! Newton-Raphson iteration : Phi_il contains current guess value
29  DO iter=1,2 ! 2 iterations should be enough
30    ! Compute pressure, A_ik
31    SELECT CASE(caldyn_thermo)
32    CASE(thermo_theta)
33      {% call() compute_p_and_Aik() %}
34        X_ij = (cpp/preff)*kappa*theta(CELL)*rho_ij
35        p_ik(CELL) = preff*(X_ij**gamma)
36        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho
37      {% endcall %}
38    CASE(thermo_entropy)
39      {% call() compute_p_and_Aik() %}
40        X_ij = log(vreff*rho_ij) + theta(CELL)/cpp
41        p_ik(CELL) = preff*exp(X_ij*gamma)
42        c2_mik = gamma*p_ik(CELL)/(rho_ij*m_ik(CELL)) ! c^2 = gamma*R*T = gamma*p/rho
43      {% endcall %}
44    CASE DEFAULT
45        PRINT *, 'caldyn_thermo not supported by compute_NH_geopot', caldyn_thermo
46        STOP
47    END SELECT
48   
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
52    SEQUENCE_C1
53      ! Compute residual R_il and B_il
54      PROLOGUE(1)
55        ! bottom interface l=1
56        ml_g2 = gm2*m_il(CELL)
57        B_il(CELL) = A_ik(CELL) + ml_g2 + tau2_g*rho_bot
58        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
59                  + tau2_g*( p_ik(CELL)-pbot+rho_bot*(Phi_il(CELL)-PHI_BOT(HIDX(CELL))) )
60      END_BLOCK
61      BODY('2,llm')
62        ! inner interfaces
63        ml_g2 = gm2*m_il(CELL)
64        B_il(CELL) = A_ik(CELL)+A_ik(DOWN(CELL)) + ml_g2
65        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
66                  + tau2_g*(p_ik(CELL)-p_ik(DOWN(CELL)))
67        ! consistency check : if Wil=0 and initial state is in hydrostatic balance
68        ! then Phi_star_il(CELL) = Phi_il(CELL) - tau^2*g^2
69        ! and residual = tau^2*(ml+(1/g)dl_pi)=0
70      END_BLOCK
71      EPILOGUE('llm+1')
72        ! top interface l=llm+1
73        ml_g2 = gm2*m_il(CELL)
74        B_il(CELL) = A_ik(DOWN(CELL)) + ml_g2
75        R_il(CELL) = ml_g2*( Phi_il(CELL)-Phi_star_il(CELL))  &
76                  + tau2_g*( ptop-p_ik(DOWN(CELL)) )
77      END_BLOCK
78      !
79      ! Forward sweep :
80      ! C(0)=0, C(l) = -A(l) / (B(l)+A(l-1)C(l-1)),
81      ! D(0)=0, D(l) = (R(l)+A(l-1)D(l-1)) / (B(l)+A(l-1)C(l-1))
82      PROLOGUE(1)
83        X_ij = 1./B_il(CELL)
84        C_ik(CELL) = -A_ik(CELL) * X_ij
85        D_il(CELL) =  R_il(CELL) * X_ij
86      END_BLOCK
87      BODY('2,llm')
88        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )
89        C_ik(CELL) = -A_ik(CELL) * X_ij
90        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij
91      END_BLOCK
92      EPILOGUE('llm+1')
93        X_ij = 1./( B_il(CELL) + A_ik(DOWN(CELL))*C_ik(DOWN(CELL)) )
94        D_il(CELL) =  (R_il(CELL)+A_ik(DOWN(CELL))*D_il(DOWN(CELL))) * X_ij
95        ! Back substitution :
96        ! x(i) = D(i)-C(i)x(i+1), x(llm+1)=0
97        ! + Newton-Raphson update
98        ! top interface l=llm+1
99        x_il(CELL) = D_il(CELL)
100        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)
101      END_BLOCK
102      BODY('llm,1,-1')
103        ! Back substitution at lower interfaces
104        x_il(CELL) = D_il(CELL) - C_ik(CELL)*x_il(UP(CELL))
105        Phi_il(CELL) = Phi_il(CELL) - x_il(CELL)
106      END_BLOCK
107    END_BLOCK
108
109    IF(debug_hevi_solver) THEN
110       PRINT *, '[hevi_solver] A,B', iter, MAXVAL(ABS(A_ik)),MAXVAL(ABS(B_il))
111       PRINT *, '[hevi_solver] C,D', iter, MAXVAL(ABS(C_ik)),MAXVAL(ABS(D_il))
112       DO l=1,llm+1
113          WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] x_il', iter,l, MAXVAL(ABS(x_il(l,:)))
114       END DO
115       DO l=1,llm+1
116          WRITE(*,'(A,I2.1,I3.2,E9.2)') '[hevi_solver] R_il', iter,l, MAXVAL(ABS(R_il(l,:)))
117       END DO
118    END IF
119
120  END DO ! Newton-Raphson
121
122  BARRIER
123
124  debug_hevi_solver=.FALSE.
125END_BLOCK
126
127KERNEL(caldyn_mil)
128  FORALL_CELLS_EXT('1', 'llm+1')
129    ON_PRIMAL
130       CST_IF(IS_BOTTOM_LEVEL,    m_il(CELL) = .5*rhodz(CELL)                     )
131       CST_IF(IS_INNER_INTERFACE, m_il(CELL) = .5*(rhodz(CELL)+rhodz(DOWN(CELL))) )
132       CST_IF(IS_TOP_INTERFACE,   m_il(CELL) = .5*rhodz(DOWN(CELL))               )
133    END_BLOCK
134  END_BLOCK
135END_BLOCK
136
137KERNEL(caldyn_solver)
138  !
139  ! Compute pressure (pres) and Exner function (pk)
140  ! kappa = R/Cp
141  ! 1-kappa = Cv/Cp
142  ! Cp/Cv = 1/(1-kappa)
143  gamma    = 1./(1.-kappa)
144  vreff    = Rd*Treff/preff ! reference specific volume
145  Cvd      = 1./(cpp-Rd)
146  Rd_preff = kappa*cpp/preff
147  FORALL_CELLS_EXT()
148    SELECT CASE(caldyn_thermo)
149    CASE(thermo_theta)
150      ON_PRIMAL
151        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL))
152        rho_ij = rho_ij*g*rhodz(CELL)
153        X_ij = Rd_preff*theta(CELL,1)*rho_ij
154        ! kappa.theta.rho = p/exner
155        ! => X = (p/p0)/(exner/Cp)
156        !      = (p/p0)^(1-kappa)
157        pres(CELL) = preff*(X_ij**gamma)  ! pressure
158        ! Compute Exner function (needed by compute_caldyn_fast)
159        ! other formulae possible if exponentiation is slow
160        pk(CELL)   = cpp*((pres(CELL)/preff)**kappa) ! Exner
161      END_BLOCK
162    CASE(thermo_entropy)
163      ON_PRIMAL
164        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL))
165        rho_ij = rho_ij*g*rhodz(CELL)
166        T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))*Cvd )
167        pres(CELL) = rho_ij*Rd*T_ij
168        pk(CELL)   = T_ij
169      END_BLOCK
170    CASE DEFAULT
171      STOP
172    END SELECT
173  END_BLOCK
174
175  ! We need a barrier here because we compute pres above and do a vertical difference below
176  BARRIER
177
178  FORALL_CELLS_EXT('1', 'llm+1')
179    ON_PRIMAL
180      CST_IFTHEN(IS_BOTTOM_LEVEL)
181        ! Lower BC
182        dW(CELL)   = (1./g)*(pbot-rho_bot*(geopot(CELL)-PHI_BOT(HIDX(CELL)))-pres(CELL)) - m_il(CELL)
183      CST_ELSEIF(IS_TOP_INTERFACE)
184        ! Top BC
185        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-ptop) - m_il(CELL)
186      CST_ELSE
187        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-pres(CELL)) - m_il(CELL)
188      CST_ENDIF
189      W(CELL)    = W(CELL)+tau*dW(CELL) ! update W
190      dPhi(CELL) = g*g*W(CELL)/m_il(CELL)
191    END_BLOCK
192  END_BLOCK
193
194  ! We need a barrier here because we update W above and do a vertical average below
195  BARRIER
196
197  FORALL_CELLS_EXT()
198    ON_PRIMAL
199      ! compute du = -0.5*g^2.grad(w^2)
200      berni(CELL) = (-.25*g*g)*((W(CELL)/m_il(CELL))**2 + (W(UP(CELL))/m_il(UP(CELL)))**2 )
201    END_BLOCK
202  END_BLOCK
203  FORALL_CELLS_EXT()
204    ON_EDGES
205      du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2))
206    END_BLOCK
207  END_BLOCK
208END_BLOCK
209
210KERNEL(caldyn_vert_NH)
211  FORALL_CELLS_EXT()
212    ON_PRIMAL
213      CST_IF(IS_BOTTOM_LEVEL,     w_ij = ( W(CELL)+.5*W(UP(CELL)) )/mass(CELL) )
214      CST_IF(IS_INNER_LAYER , w_ij = .5*( W(CELL)+W(UP(CELL)) )/mass(CELL) )
215      CST_IF(IS_TOP_LAYER,        w_ij = ( .5*W(CELL)+W(UP(CELL)) )/mass(CELL) )
216      wflux_ij = .5*(wflux(CELL)+wflux(UP(CELL)))
217      W_etadot(CELL) = wflux_ij*w_ij
218      eta_dot(CELL) = wflux_ij / mass(CELL)
219      wcov(CELL) = w_ij*(geopot(UP(CELL))-geopot(CELL))
220    END_BLOCK
221  END_BLOCK
222  ! add NH term to du
223  FORALL_CELLS()
224    ON_EDGES
225      du(EDGE) = du(EDGE) - .5*(wcov(CELL1)+wcov(CELL2))*SIGN*(eta_dot(CELL2)-eta_dot(CELL1))
226    END_BLOCK
227  END_BLOCK
228  ! add NH terms to dW, dPhi
229  ! FIXME : TODO top and bottom
230  FORALL_CELLS('2','llm')
231    ON_PRIMAL
232      dPhi(CELL)=dPhi(CELL)-wflux(CELL)*(geopot(UP(CELL))-geopot(DOWN(CELL)))/(mass(DOWN(CELL))+mass(CELL))
233    END_BLOCK
234  END_BLOCK
235
236  ! We need a barrier here because we compute W_etadot above and do a vertical difference below
237  BARRIER
238
239  FORALL_CELLS('1','llm+1')
240    ON_PRIMAL
241      CST_IFTHEN(IS_BOTTOM_LEVEL)
242         dW(CELL) = dW(CELL) - W_etadot(CELL)
243      CST_ELSEIF(IS_TOP_INTERFACE)
244         dW(CELL) = dW(CELL) + W_etadot(DOWN(CELL))
245      CST_ELSE
246         dW(CELL) = dW(CELL) + W_etadot(DOWN(CELL)) - W_etadot(CELL)
247      CST_ENDIF
248    END_BLOCK
249  END_BLOCK
250END_BLOCK
251
Note: See TracBrowser for help on using the repository browser.