source: codes/icosagcm/devel/src/dynamics/compute_caldyn_solver.F90

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

devel : reduced ground stiffness for better stability (NH)
File size: 7.9 KB
Line 
1MODULE compute_caldyn_solver_mod
2  USE grid_param
3  USE compute_NH_geopot_mod, ONLY : compute_NH_geopot, pbot, rho_bot
4  IMPLICIT NONE
5  PRIVATE
6
7#include "../unstructured/unstructured.h90"
8
9  PUBLIC :: compute_caldyn_solver
10
11CONTAINS
12
13#ifdef BEGIN_DYSL
14
15KERNEL(caldyn_solver)
16  !   
17  ! Compute pressure (pres) and Exner function (pk)
18  ! kappa = R/Cp
19  ! 1-kappa = Cv/Cp
20  ! Cp/Cv = 1/(1-kappa)
21  gamma    = 1./(1.-kappa)
22  vreff    = Rd*Treff/preff ! reference specific volume
23  Cvd      = 1./(cpp-Rd)
24  Rd_preff = kappa*cpp/preff
25  FORALL_CELLS_EXT()
26    SELECT CASE(caldyn_thermo)
27    CASE(thermo_theta)
28      ON_PRIMAL
29        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL))
30        rho_ij = rho_ij*g*rhodz(CELL)
31        X_ij = Rd_preff*theta(CELL,1)*rho_ij
32        ! kappa.theta.rho = p/exner
33        ! => X = (p/p0)/(exner/Cp)
34        !      = (p/p0)^(1-kappa)
35        pres(CELL) = preff*(X_ij**gamma)  ! pressure
36        ! Compute Exner function (needed by compute_caldyn_fast)
37        ! other formulae possible if exponentiation is slow
38        pk(CELL)   = cpp*((pres(CELL)/preff)**kappa) ! Exner
39      END_BLOCK
40    CASE(thermo_entropy)
41      ON_PRIMAL
42        rho_ij = 1./(geopot(UP(CELL))-geopot(CELL))
43        rho_ij = rho_ij*g*rhodz(CELL)
44        T_ij = Treff*exp( (theta(CELL,1)+Rd*log(vreff*rho_ij))*Cvd )
45        pres(CELL) = rho_ij*Rd*T_ij
46        pk(CELL)   = T_ij
47      END_BLOCK
48    CASE DEFAULT
49      STOP
50    END SELECT
51  END_BLOCK
52
53  ! We need a barrier here because we compute pres above and do a vertical difference below
54  BARRIER
55
56  FORALL_CELLS_EXT('1', 'llm+1')
57    ON_PRIMAL
58      CST_IFTHEN(IS_BOTTOM_LEVEL)
59        ! Lower BC
60        dW(CELL)   = (1./g)*(pbot-rho_bot*(geopot(CELL)-PHI_BOT(HIDX(CELL)))-pres(CELL)) - m_il(CELL)
61      CST_ELSEIF(IS_TOP_INTERFACE)
62        ! Top BC
63        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-ptop) - m_il(CELL)
64      CST_ELSE
65        dW(CELL)   = (1./g)*(pres(DOWN(CELL))-pres(CELL)) - m_il(CELL)
66      CST_ENDIF
67      W(CELL)    = W(CELL)+tau*dW(CELL) ! update W
68      dPhi(CELL) = g*g*W(CELL)/m_il(CELL)
69    END_BLOCK
70  END_BLOCK
71
72  ! We need a barrier here because we update W above and do a vertical average below
73  BARRIER
74
75  FORALL_CELLS_EXT()
76    ON_PRIMAL
77      ! compute du = -0.5*g^2.grad(w^2)
78      berni(CELL) = (-.25*g*g)*((W(CELL)/m_il(CELL))**2 + (W(UP(CELL))/m_il(UP(CELL)))**2 )
79    END_BLOCK
80  END_BLOCK
81  FORALL_CELLS_EXT()
82    ON_EDGES
83      du(EDGE) = SIGN*(berni(CELL1)-berni(CELL2))
84    END_BLOCK
85  END_BLOCK
86END_BLOCK
87
88#endif END_DYSL
89
90SUBROUTINE compute_caldyn_solver_unst(tau,rhodz,theta, berni,pres,m_il, pk,geopot,W,dPhi,dW,du)
91  USE ISO_C_BINDING, only : C_DOUBLE, C_FLOAT
92  USE earth_const
93  USE trace
94  USE grid_param, ONLY : nqdyn
95  USE disvert_mod, ONLY : ptop
96  USE data_unstructured_mod, ONLY : enter_trace, exit_trace, &
97       id_solver, left, right, PHI_BOT
98   
99  USE compute_NH_geopot_mod, ONLY : compute_NH_geopot_unst
100  NUM, INTENT(IN) :: tau 
101  FIELD_MASS   :: rhodz,pk,berni,pres    ! IN, OUT, BUF*2
102  FIELD_THETA  :: theta                  ! IN
103  FIELD_GEOPOT :: geopot,W,dPhi,dW, m_il ! INOUT,INOUT, OUT,OUT, BUF
104  FIELD_U      :: du                     ! OUT
105  DECLARE_INDICES
106  NUM :: X_ij, rho_ij, T_ij, gamma, Cvd, vreff, Rd_preff
107#define PHI_BOT(ij) Phi_bot
108#include "../kernels_unst/caldyn_mil.k90"
109  IF(tau>0) THEN ! solve implicit problem for geopotential
110    CALL compute_NH_geopot_unst(tau, rhodz, m_il, theta, W, geopot)
111  END IF
112  START_TRACE(id_solver, 7,0,1)
113#include "../kernels_unst/caldyn_solver.k90"
114  STOP_TRACE
115#undef PHI_BOT
116END SUBROUTINE compute_caldyn_solver_unst
117
118  SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)
119    USE icosa
120    USE caldyn_vars_mod
121    USE trace
122    USE omp_para, ONLY : ll_begin, ll_end,ll_beginp1,ll_endp1
123    USE disvert_mod, ONLY : ptop
124    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
125    REAL(rstd),INTENT(IN)    :: phis(iim*jjm)
126    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
127    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
128    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)
129    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
130    REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0
131    REAL(rstd),INTENT(OUT)   :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces
132    REAL(rstd),INTENT(OUT)   :: pres(iim*jjm,llm)          ! pressure
133    REAL(rstd),INTENT(OUT)   :: dW(iim*jjm,llm+1)
134    REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1)
135    REAL(rstd),INTENT(OUT)   :: du(3*iim*jjm,llm)
136
137    REAL(rstd) :: berni(iim*jjm,llm)         ! (W/m_il)^2
138    REAL(rstd) :: berni1(iim*jjm)         ! (W/m_il)^2
139    !REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd, Rd_preff
140    REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Cvd, Rd_preff
141    INTEGER    :: ij, l
142
143    CALL trace_start("compute_caldyn_solver")
144
145    !Rd=cpp*kappa
146
147    IF(dysl) THEN
148
149!$OMP BARRIER
150
151#include "../kernels_hex/caldyn_mil.k90"
152  IF(tau>0) THEN ! solve implicit problem for geopotential
153    CALL compute_NH_geopot(tau,phis, rhodz, m_il, theta, W, geopot)
154  END IF
155#define PHI_BOT(ij) phis(ij)
156#include "../kernels_hex/caldyn_solver.k90"
157#undef PHI_BOT
158!$OMP BARRIER
159
160    ELSE
161
162#define BERNI(ij) berni1(ij)
163    ! FIXME : vertical OpenMP parallelism will not work
164
165    ! average m_ik to interfaces => m_il
166    !DIR$ SIMD
167    DO ij=ij_begin_ext,ij_end_ext
168       m_il(ij,1) = .5*rhodz(ij,1)
169    ENDDO
170    DO l=2,llm
171       !DIR$ SIMD
172       DO ij=ij_begin_ext,ij_end_ext
173          m_il(ij,l) = .5*(rhodz(ij,l-1)+rhodz(ij,l))
174       ENDDO
175    ENDDO
176    !DIR$ SIMD
177    DO ij=ij_begin_ext,ij_end_ext
178       m_il(ij,llm+1) = .5*rhodz(ij,llm)
179    ENDDO
180
181    IF(tau>0) THEN ! solve implicit problem for geopotential
182       CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot)
183    END IF
184
185    ! Compute pressure, stored temporarily in pk
186    ! kappa = R/Cp
187    ! 1-kappa = Cv/Cp
188    ! Cp/Cv = 1/(1-kappa)
189    gamma = 1./(1.-kappa)
190    DO l=1,llm
191       !DIR$ SIMD
192       DO ij=ij_begin_ext,ij_end_ext
193          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))
194          X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij
195          ! kappa.theta.rho = p/exner
196          ! => X = (p/p0)/(exner/Cp)
197          !      = (p/p0)^(1-kappa)
198          pk(ij,l) = preff*(X_ij**gamma)
199       ENDDO
200    ENDDO
201
202    ! Update W, compute tendencies
203    DO l=2,llm
204       !DIR$ SIMD
205       DO ij=ij_begin_ext,ij_end_ext
206          dW(ij,l)   = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il(ij,l)
207          W(ij,l)    = W(ij,l)+tau*dW(ij,l) ! update W
208          dPhi(ij,l) = g*g*W(ij,l)/m_il(ij,l)
209       ENDDO
210       !       PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l)))
211       !       PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l)))
212    ENDDO
213    ! Lower BC (FIXME : no orography yet !)
214    DO ij=ij_begin,ij_end         
215       dPhi(ij,1)=0
216       W(ij,1)=0
217       dW(ij,1)=0
218       dPhi(ij,llm+1)=0 ! rigid lid
219       W(ij,llm+1)=0
220       dW(ij,llm+1)=0
221    ENDDO
222    ! Upper BC p=ptop
223    !    DO ij=ij_omp_begin_ext,ij_omp_end_ext
224    !       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm)
225    !       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm)
226    !    ENDDO
227   
228    ! Compute Exner function (needed by compute_caldyn_fast) and du=-g^2.grad(w^2)
229    DO l=1,llm
230       !DIR$ SIMD
231       DO ij=ij_begin_ext,ij_end_ext
232          pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow
233          BERNI(ij) = (-.25*g*g)*(              &
234                 (W(ij,l)/m_il(ij,l))**2       &
235               + (W(ij,l+1)/m_il(ij,l+1))**2 )
236       ENDDO
237       DO ij=ij_begin,ij_end         
238          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
239          du(ij+u_lup,l)   = ne_lup  *(BERNI(ij)-BERNI(ij+t_lup))
240          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
241       ENDDO
242    ENDDO
243#undef BERNI
244
245    END IF ! dysl
246
247    CALL trace_end("compute_caldyn_solver")
248   
249  END SUBROUTINE compute_caldyn_solver
250 
251END MODULE compute_caldyn_solver_mod
Note: See TracBrowser for help on using the repository browser.