Changeset 538 for codes/icosagcm/devel
- Timestamp:
- 06/09/17 16:13:47 (7 years ago)
- Location:
- codes/icosagcm/devel/src/dynamics
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/dynamics/caldyn_kernels_base.F90
r537 r538 49 49 50 50 #ifdef CPP_DYSL 51 !#if 0 51 52 #include "../kernels/compute_geopot.k90" 52 53 #else … … 303 304 REAL(rstd),INTENT(INOUT) :: dW(iim*jjm,llm+1) 304 305 ! local arrays 305 REAL(rstd) :: eta_dot(iim*jjm ) ! eta_dot in full layers306 REAL(rstd) :: wcov(iim*jjm ) ! covariant vertical momentum in full layers306 REAL(rstd) :: eta_dot(iim*jjm, llm) ! eta_dot in full layers 307 REAL(rstd) :: wcov(iim*jjm,llm) ! covariant vertical momentum in full layers 307 308 REAL(rstd) :: W_etadot(iim*jjm,llm) ! vertical flux of vertical momentum 308 309 ! indices and temporary values … … 312 313 CALL trace_start("compute_caldyn_vert_nh") 313 314 315 #ifdef CPP_DYSL 316 !#if 0 317 #include "../kernels/caldyn_vert_NH.k90" 318 #else 319 #define ETA_DOT(ij) eta_dot(ij,1) 320 #define WCOV(ij) wcov(ij,1) 321 314 322 DO l=ll_begin,ll_end 315 323 ! compute the local arrays … … 319 327 w_ij = .5*(W(ij,l)+W(ij,l+1))/mass(ij,l) 320 328 W_etadot(ij,l) = wflux_ij*w_ij 321 eta_dot(ij) = wflux_ij / mass(ij,l)322 wcov(ij) = w_ij*(geopot(ij,l+1)-geopot(ij,l))329 ETA_DOT(ij) = wflux_ij / mass(ij,l) 330 WCOV(ij) = w_ij*(geopot(ij,l+1)-geopot(ij,l)) 323 331 ENDDO 324 332 ! add NH term to du … … 326 334 DO ij=ij_begin,ij_end 327 335 du(ij+u_right,l) = du(ij+u_right,l) & 328 - .5*( wcov(ij+t_right)+wcov(ij)) &329 *ne_right*( eta_dot(ij+t_right)-eta_dot(ij))336 - .5*(WCOV(ij+t_right)+WCOV(ij)) & 337 *ne_right*(ETA_DOT(ij+t_right)-ETA_DOT(ij)) 330 338 du(ij+u_lup,l) = du(ij+u_lup,l) & 331 - .5*( wcov(ij+t_lup)+wcov(ij)) &332 *ne_lup*( eta_dot(ij+t_lup)-eta_dot(ij))339 - .5*(WCOV(ij+t_lup)+WCOV(ij)) & 340 *ne_lup*(ETA_DOT(ij+t_lup)-ETA_DOT(ij)) 333 341 du(ij+u_ldown,l) = du(ij+u_ldown,l) & 334 - .5*( wcov(ij+t_ldown)+wcov(ij)) &335 *ne_ldown*( eta_dot(ij+t_ldown)-eta_dot(ij))342 - .5*(WCOV(ij+t_ldown)+WCOV(ij)) & 343 *ne_ldown*(ETA_DOT(ij+t_ldown)-ETA_DOT(ij)) 336 344 END DO 337 345 ENDDO … … 352 360 END DO 353 361 END DO 362 363 #undef ETA_DOT 364 #undef WCOV 365 #endif 366 354 367 CALL trace_end("compute_caldyn_vert_nh") 355 368 -
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.