Changeset 853 for codes/icosagcm/devel
- Timestamp:
- 05/05/19 16:17:23 (5 years ago)
- Location:
- codes/icosagcm/devel/src/dynamics
- Files:
-
- 1 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90
r850 r853 10 10 USE compute_caldyn_slow_hydro_mod, ONLY : compute_caldyn_slow_hydro 11 11 USE compute_caldyn_slow_NH_mod, ONLY : compute_caldyn_slow_NH 12 USE compute_caldyn_solver_mod, ONLY : compute_caldyn_solver 12 13 IMPLICIT NONE 13 14 PRIVATE -
codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
r850 r853 13 13 LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 14 14 15 PUBLIC :: compute_caldyn_ solver, compute_caldyn_fast15 PUBLIC :: compute_caldyn_fast,compute_NH_geopot 16 16 17 17 CONTAINS … … 183 183 END SUBROUTINE compute_NH_geopot 184 184 185 SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)186 REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs187 REAL(rstd),INTENT(IN) :: phis(iim*jjm)188 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm)189 REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn)190 REAL(rstd),INTENT(OUT) :: pk(iim*jjm,llm)191 REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)192 REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0193 REAL(rstd),INTENT(OUT) :: m_il(iim*jjm,llm+1) ! rhodz averaged to interfaces194 REAL(rstd),INTENT(OUT) :: pres(iim*jjm,llm) ! pressure195 REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1)196 REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1)197 REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)198 199 REAL(rstd) :: berni(iim*jjm,llm) ! (W/m_il)^2200 REAL(rstd) :: berni1(iim*jjm) ! (W/m_il)^2201 REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd, Rd_preff202 INTEGER :: ij, l203 204 CALL trace_start("compute_caldyn_solver")205 206 Rd=cpp*kappa207 208 IF(dysl) THEN209 210 !$OMP BARRIER211 212 #include "../kernels_hex/caldyn_mil.k90"213 IF(tau>0) THEN ! solve implicit problem for geopotential214 CALL compute_NH_geopot(tau,phis, rhodz, m_il, theta, W, geopot)215 END IF216 #define PHI_BOT(ij) phis(ij)217 #include "../kernels_hex/caldyn_solver.k90"218 #undef PHI_BOT219 !$OMP BARRIER220 221 ELSE222 223 #define BERNI(ij) berni1(ij)224 ! FIXME : vertical OpenMP parallelism will not work225 226 ! average m_ik to interfaces => m_il227 !DIR$ SIMD228 DO ij=ij_begin_ext,ij_end_ext229 m_il(ij,1) = .5*rhodz(ij,1)230 ENDDO231 DO l=2,llm232 !DIR$ SIMD233 DO ij=ij_begin_ext,ij_end_ext234 m_il(ij,l) = .5*(rhodz(ij,l-1)+rhodz(ij,l))235 ENDDO236 ENDDO237 !DIR$ SIMD238 DO ij=ij_begin_ext,ij_end_ext239 m_il(ij,llm+1) = .5*rhodz(ij,llm)240 ENDDO241 242 IF(tau>0) THEN ! solve implicit problem for geopotential243 CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot)244 END IF245 246 ! Compute pressure, stored temporarily in pk247 ! kappa = R/Cp248 ! 1-kappa = Cv/Cp249 ! Cp/Cv = 1/(1-kappa)250 gamma = 1./(1.-kappa)251 DO l=1,llm252 !DIR$ SIMD253 DO ij=ij_begin_ext,ij_end_ext254 rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))255 X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij256 ! kappa.theta.rho = p/exner257 ! => X = (p/p0)/(exner/Cp)258 ! = (p/p0)^(1-kappa)259 pk(ij,l) = preff*(X_ij**gamma)260 ENDDO261 ENDDO262 263 ! Update W, compute tendencies264 DO l=2,llm265 !DIR$ SIMD266 DO ij=ij_begin_ext,ij_end_ext267 dW(ij,l) = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il(ij,l)268 W(ij,l) = W(ij,l)+tau*dW(ij,l) ! update W269 dPhi(ij,l) = g*g*W(ij,l)/m_il(ij,l)270 ENDDO271 ! PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l)))272 ! PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l)))273 ENDDO274 ! Lower BC (FIXME : no orography yet !)275 DO ij=ij_begin,ij_end276 dPhi(ij,1)=0277 W(ij,1)=0278 dW(ij,1)=0279 dPhi(ij,llm+1)=0 ! rigid lid280 W(ij,llm+1)=0281 dW(ij,llm+1)=0282 ENDDO283 ! Upper BC p=ptop284 ! DO ij=ij_omp_begin_ext,ij_omp_end_ext285 ! dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm)286 ! dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm)287 ! ENDDO288 289 ! Compute Exner function (needed by compute_caldyn_fast) and du=-g^2.grad(w^2)290 DO l=1,llm291 !DIR$ SIMD292 DO ij=ij_begin_ext,ij_end_ext293 pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow294 BERNI(ij) = (-.25*g*g)*( &295 (W(ij,l)/m_il(ij,l))**2 &296 + (W(ij,l+1)/m_il(ij,l+1))**2 )297 ENDDO298 DO ij=ij_begin,ij_end299 du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))300 du(ij+u_lup,l) = ne_lup *(BERNI(ij)-BERNI(ij+t_lup))301 du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))302 ENDDO303 ENDDO304 #undef BERNI305 306 END IF ! dysl307 308 CALL trace_end("compute_caldyn_solver")309 310 END SUBROUTINE compute_caldyn_solver311 312 185 SUBROUTINE compute_caldyn_fast(tau,u,rhodz,theta,pk,geopot,du) 313 186 REAL(rstd),INTENT(IN) :: tau ! "solve" u-tau*du/dt = rhs
Note: See TracChangeset
for help on using the changeset viewer.