Changeset 538 for codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
- Timestamp:
- 06/09/17 16:13:47 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
r537 r538 11 11 REAL(rstd), PARAMETER :: pbot=1e5, Phi_bot=0., rho_bot=1e6 ! FIXME 12 12 13 LOGICAL, PARAMETER :: debug_hevi_solver = .FALSE. 13 LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 14 LOGICAL, PARAMETER :: rigid=.TRUE. 14 15 15 16 PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, & … … 117 118 REAL(rstd) :: wil, tau2_g, g2, gm2, ml_g2, c2_mik 118 119 119 INTEGER :: iter, ij, l 120 INTEGER :: iter, ij, l, ij_omp_begin_ext, ij_omp_end_ext 121 122 CALL distrib_level(ij_end_ext-ij_begin_ext+1,ij_omp_begin_ext,ij_omp_end_ext) 123 ij_omp_begin_ext=ij_omp_begin_ext+ij_begin_ext-1 124 ij_omp_end_ext=ij_omp_end_ext+ij_begin_ext-1 125 126 #ifdef CPP_DYSL 127 !#if 0 128 #include "../kernels/compute_NH_geopot.k90" 129 #else 120 130 121 131 ! FIXME : vertical OpenMP parallelism will not work … … 251 261 252 262 END DO ! Newton-Raphson 263 264 #endif 253 265 254 266 END SUBROUTINE compute_NH_geopot … … 257 269 REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs 258 270 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) 259 REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm )271 REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm,nqdyn) 260 272 REAL(rstd),INTENT(OUT) :: pk(iim*jjm,llm) 261 273 REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1) … … 266 278 267 279 REAL(rstd) :: m_il(iim*jjm,llm+1) ! rhodz averaged to interfaces 268 REAL(rstd) :: berni(iim*jjm) ! (W/m_il)^2 269 REAL(rstd) :: gamma, rho_ij, X_ij, Y_ij 280 REAL(rstd) :: pres(iim*jjm,llm) ! pressure 281 REAL(rstd) :: berni(iim*jjm,llm) ! (W/m_il)^2 282 REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd 270 283 INTEGER :: ij, l 271 284 272 285 CALL trace_start("compute_caldyn_solver") 273 286 287 Rd=cpp*kappa 288 289 #ifdef CPP_DYSL 290 !#if 0 291 #include "../kernels/caldyn_solver.k90" 292 #else 293 #define BERNI(ij) berni(ij,1) 274 294 ! FIXME : vertical OpenMP parallelism will not work 275 295 … … 303 323 DO ij=ij_begin_ext,ij_end_ext 304 324 rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l)) 305 X_ij = (cpp/preff)*kappa*theta(ij,l )*rho_ij325 X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij 306 326 ! kappa.theta.rho = p/exner 307 327 ! => X = (p/p0)/(exner/Cp) … … 342 362 DO ij=ij_begin_ext,ij_end_ext 343 363 pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow 344 berni(ij) = (-.25*g*g)*( &364 BERNI(ij) = (-.25*g*g)*( & 345 365 (W(ij,l)/m_il(ij,l))**2 & 346 366 + (W(ij,l+1)/m_il(ij,l+1))**2 ) 347 367 ENDDO 348 368 DO ij=ij_begin,ij_end 349 du(ij+u_right,l) = ne_right*( berni(ij)-berni(ij+t_right))350 du(ij+u_lup,l) = ne_lup *( berni(ij)-berni(ij+t_lup))351 du(ij+u_ldown,l) = ne_ldown*( berni(ij)-berni(ij+t_ldown))369 du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right)) 370 du(ij+u_lup,l) = ne_lup *(BERNI(ij)-BERNI(ij+t_lup)) 371 du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown)) 352 372 ENDDO 353 373 ENDDO 354 374 #undef BERNI 375 #endif 376 355 377 CALL trace_end("compute_caldyn_solver") 356 378 … … 488 510 REAL(rstd),INTENT(INOUT) :: du(3*iim*jjm,llm) 489 511 490 REAL(rstd) :: Ftheta(3*iim*jjm ) ! potential temperature flux491 REAL(rstd) :: uu_right, uu_lup, uu_ldown 512 REAL(rstd) :: Ftheta(3*iim*jjm,llm) ! potential temperature flux 513 REAL(rstd) :: uu_right, uu_lup, uu_ldown, du_trisk, divF 492 514 INTEGER :: ij,iq,l,kdown 493 515 494 516 CALL trace_start("compute_caldyn_Coriolis") 517 518 #ifdef CPP_DYSL 519 !#if 0 520 #include "../kernels/coriolis.k90" 521 #else 522 #define FTHETA(ij) Ftheta(ij,1) 495 523 496 524 DO l=ll_begin, ll_end … … 499 527 !DIR$ SIMD 500 528 DO ij=ij_begin_ext,ij_end_ext 501 F theta(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) &529 FTHETA(ij+u_right) = 0.5*(theta(ij,l,iq)+theta(ij+t_right,l,iq)) & 502 530 * hflux(ij+u_right,l) 503 F theta(ij+u_lup) = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) &531 FTHETA(ij+u_lup) = 0.5*(theta(ij,l,iq)+theta(ij+t_lup,l,iq)) & 504 532 * hflux(ij+u_lup,l) 505 F theta(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) &533 FTHETA(ij+u_ldown) = 0.5*(theta(ij,l,iq)+theta(ij+t_ldown,l,iq)) & 506 534 * hflux(ij+u_ldown,l) 507 535 END DO … … 511 539 ! dtheta_rhodz = -div(flux.theta) 512 540 dtheta_rhodz(ij,l,iq)= & 513 -1./Ai(ij)*(ne_right*F theta(ij+u_right) + &514 ne_rup*F theta(ij+u_rup) + &515 ne_lup*F theta(ij+u_lup) + &516 ne_left*F theta(ij+u_left) + &517 ne_ldown*F theta(ij+u_ldown) + &518 ne_rdown*F theta(ij+u_rdown) )541 -1./Ai(ij)*(ne_right*FTHETA(ij+u_right) + & 542 ne_rup*FTHETA(ij+u_rup) + & 543 ne_lup*FTHETA(ij+u_lup) + & 544 ne_left*FTHETA(ij+u_left) + & 545 ne_ldown*FTHETA(ij+u_ldown) + & 546 ne_rdown*FTHETA(ij+u_rdown) ) 519 547 END DO 520 548 END DO … … 626 654 STOP 627 655 END SELECT 656 #undef FTHETA 657 #endif 628 658 629 659 CALL trace_end("compute_caldyn_Coriolis") … … 645 675 646 676 #ifdef CPP_DYSL 677 !#if 0 647 678 #define BERNI(ij,l) berni(ij,l) 648 679 #include "../kernels/caldyn_slow_hydro.k90" 680 #undef BERNI 649 681 #else 650 682 #define BERNI(ij) berni(ij,1) … … 696 728 END IF 697 729 END DO 730 #undef BERNI 698 731 #endif 699 #undef BERNI700 732 CALL trace_end("compute_caldyn_slow_hydro") 701 733 END SUBROUTINE compute_caldyn_slow_hydro
Note: See TracChangeset
for help on using the changeset viewer.