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

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

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