Changeset 850
- Timestamp:
- 05/05/19 14:10:43 (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
r849 r850 9 9 USE compute_caldyn_Coriolis_mod, ONLY : compute_caldyn_Coriolis 10 10 USE compute_caldyn_slow_hydro_mod, ONLY : compute_caldyn_slow_hydro 11 USE compute_caldyn_slow_NH_mod, ONLY : compute_caldyn_slow_NH 11 12 IMPLICIT NONE 12 13 PRIVATE -
codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90
r849 r850 13 13 LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 14 14 15 PUBLIC :: compute_caldyn_slow_NH, & 16 compute_caldyn_solver, compute_caldyn_fast 15 PUBLIC :: compute_caldyn_solver, compute_caldyn_fast 17 16 18 17 CONTAINS … … 434 433 END SUBROUTINE compute_caldyn_fast 435 434 436 437 SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)438 REAL(rstd),INTENT(IN) :: u(3*iim*jjm,llm) ! prognostic "velocity"439 REAL(rstd),INTENT(IN) :: rhodz(iim*jjm,llm) ! rho*dz440 REAL(rstd),INTENT(IN) :: Phi(iim*jjm,llm+1) ! prognostic geopotential441 REAL(rstd),INTENT(IN) :: W(iim*jjm,llm+1) ! prognostic vertical momentum442 443 REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s444 REAL(rstd),INTENT(OUT) :: du(3*iim*jjm,llm)445 REAL(rstd),INTENT(OUT) :: dW(iim*jjm,llm+1)446 REAL(rstd),INTENT(OUT) :: dPhi(iim*jjm,llm+1)447 448 REAL(rstd) :: w_il(iim*jjm,llm+1) ! Wil/mil449 REAL(rstd) :: F_el(3*iim*jjm,llm+1) ! NH mass flux450 REAL(rstd) :: gradPhi2(iim*jjm,llm+1) ! grad_Phi**2451 REAL(rstd) :: DePhil(3*iim*jjm,llm+1) ! grad(Phi)452 453 INTEGER :: ij,l,kdown,kup454 REAL(rstd) :: W_el, W2_el, uu_right, uu_lup, uu_ldown, gPhi2, dP, divG, u2, uu455 456 REAL(rstd) :: berni(iim*jjm,llm) ! Bernoulli function457 REAL(rstd) :: G_el(3*iim*jjm,llm+1) ! horizontal flux of W458 REAL(rstd) :: v_el(3*iim*jjm,llm+1)459 460 REAL(rstd) :: berni1(iim*jjm) ! Bernoulli function461 REAL(rstd) :: G_el1(3*iim*jjm) ! horizontal flux of W462 REAL(rstd) :: v_el1(3*iim*jjm)463 464 CALL trace_start("compute_caldyn_slow_NH")465 466 IF(dysl) THEN467 468 !$OMP BARRIER469 #include "../kernels_hex/caldyn_slow_NH.k90"470 !$OMP BARRIER471 472 ELSE473 474 #define BERNI(ij) berni1(ij)475 #define G_EL(ij) G_el1(ij)476 #define V_EL(ij) v_el1(ij)477 478 DO l=ll_begin, ll_endp1 ! compute on l levels (interfaces)479 IF(l==1) THEN480 kdown=1481 ELSE482 kdown=l-1483 END IF484 IF(l==llm+1) THEN485 kup=llm486 ELSE487 kup=l488 END IF489 ! below : "checked" means "formula also valid when kup=kdown (top/bottom)"490 ! compute mil, wil=Wil/mil491 DO ij=ij_begin_ext, ij_end_ext492 w_il(ij,l) = 2.*W(ij,l)/(rhodz(ij,kdown)+rhodz(ij,kup)) ! checked493 END DO494 ! compute DePhi, v_el, G_el, F_el495 ! v_el, W2_el and therefore G_el incorporate metric factor le_de496 ! while DePhil, W_el and F_el don't497 DO ij=ij_begin_ext, ij_end_ext498 ! Compute on edge 'right'499 W_el = .5*( W(ij,l)+W(ij+t_right,l) )500 DePhil(ij+u_right,l) = ne_right*(Phi(ij+t_right,l)-Phi(ij,l))501 F_el(ij+u_right,l) = DePhil(ij+u_right,l)*W_el502 W2_el = .5*le_de(ij+u_right) * &503 ( W(ij,l)*w_il(ij,l) + W(ij+t_right,l)*w_il(ij+t_right,l) )504 V_EL(ij+u_right) = .5*le_de(ij+u_right)*(u(ij+u_right,kup)+u(ij+u_right,kdown)) ! checked505 G_EL(ij+u_right) = V_EL(ij+u_right)*W_el - DePhil(ij+u_right,l)*W2_el506 ! Compute on edge 'lup'507 W_el = .5*( W(ij,l)+W(ij+t_lup,l) )508 DePhil(ij+u_lup,l) = ne_lup*(Phi(ij+t_lup,l)-Phi(ij,l))509 F_el(ij+u_lup,l) = DePhil(ij+u_lup,l)*W_el510 W2_el = .5*le_de(ij+u_lup) * &511 ( W(ij,l)*w_il(ij,l) + W(ij+t_lup,l)*w_il(ij+t_lup,l) )512 V_EL(ij+u_lup) = .5*le_de(ij+u_lup)*( u(ij+u_lup,kup) + u(ij+u_lup,kdown)) ! checked513 G_EL(ij+u_lup) = V_EL(ij+u_lup)*W_el - DePhil(ij+u_lup,l)*W2_el514 ! Compute on edge 'ldown'515 W_el = .5*( W(ij,l)+W(ij+t_ldown,l) )516 DePhil(ij+u_ldown,l) = ne_ldown*(Phi(ij+t_ldown,l)-Phi(ij,l))517 F_el(ij+u_ldown,l) = DePhil(ij+u_ldown,l)*W_el518 W2_el = .5*le_de(ij+u_ldown) * &519 ( W(ij,l)*w_il(ij,l) + W(ij+t_ldown,l)*w_il(ij+t_ldown,l) )520 V_EL(ij+u_ldown) = .5*le_de(ij+u_ldown)*( u(ij+u_ldown,kup) + u(ij+u_ldown,kdown)) ! checked521 G_EL(ij+u_ldown) = V_EL(ij+u_ldown)*W_el - DePhil(ij+u_ldown,l)*W2_el522 END DO523 ! compute GradPhi2, dPhi, dW524 DO ij=ij_begin_ext, ij_end_ext525 gradPhi2(ij,l) = &526 1/(2*Ai(ij))*(le_de(ij+u_right)*DePhil(ij+u_right,l)**2 + &527 le_de(ij+u_rup)*DePhil(ij+u_rup,l)**2 + &528 le_de(ij+u_lup)*DePhil(ij+u_lup,l)**2 + &529 le_de(ij+u_left)*DePhil(ij+u_left,l)**2 + &530 le_de(ij+u_ldown)*DePhil(ij+u_ldown,l)**2 + &531 le_de(ij+u_rdown)*DePhil(ij+u_rdown,l)**2 )532 533 dPhi(ij,l) = gradPhi2(ij,l)*w_il(ij,l) -1/(2*Ai(ij))* &534 ( DePhil(ij+u_right,l)*V_EL(ij+u_right) + & ! -v.gradPhi,535 DePhil(ij+u_rup,l)*V_EL(ij+u_rup) + & ! v_el already has le_de536 DePhil(ij+u_lup,l)*V_EL(ij+u_lup) + &537 DePhil(ij+u_left,l)*V_EL(ij+u_left) + &538 DePhil(ij+u_ldown,l)*V_EL(ij+u_ldown) + &539 DePhil(ij+u_rdown,l)*V_EL(ij+u_rdown) )540 541 dW(ij,l) = -1./Ai(ij)*( & ! -div(G_el),542 ne_right*G_EL(ij+u_right) + & ! G_el already has le_de543 ne_rup*G_EL(ij+u_rup) + &544 ne_lup*G_EL(ij+u_lup) + &545 ne_left*G_EL(ij+u_left) + &546 ne_ldown*G_EL(ij+u_ldown) + &547 ne_rdown*G_EL(ij+u_rdown))548 END DO549 END DO550 551 DO l=ll_begin, ll_end ! compute on k levels (layers)552 ! Compute berni at scalar points553 DO ij=ij_begin_ext, ij_end_ext554 BERNI(ij) = &555 1/(4*Ai(ij))*( &556 le_de(ij+u_right)*u(ij+u_right,l)**2 + &557 le_de(ij+u_rup)*u(ij+u_rup,l)**2 + &558 le_de(ij+u_lup)*u(ij+u_lup,l)**2 + &559 le_de(ij+u_left)*u(ij+u_left,l)**2 + &560 le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 + &561 le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 ) &562 - .25*( gradPhi2(ij,l) *w_il(ij,l)**2 + &563 gradPhi2(ij,l+1)*w_il(ij,l+1)**2 )564 END DO565 ! Compute mass flux and grad(berni) at edges566 DO ij=ij_begin_ext, ij_end_ext567 ! Compute on edge 'right'568 uu_right = 0.5*(rhodz(ij,l)+rhodz(ij+t_right,l))*u(ij+u_right,l) &569 -0.5*(F_el(ij+u_right,l)+F_el(ij+u_right,l+1))570 hflux(ij+u_right,l) = uu_right*le_de(ij+u_right)571 du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))572 ! Compute on edge 'lup'573 uu_lup = 0.5*(rhodz(ij,l)+rhodz(ij+t_lup,l))*u(ij+u_lup,l) &574 -0.5*(F_el(ij+u_lup,l)+F_el(ij+u_lup,l+1))575 hflux(ij+u_lup,l) = uu_lup*le_de(ij+u_lup)576 du(ij+u_lup,l) = ne_lup*(BERNI(ij)-BERNI(ij+t_lup))577 ! Compute on edge 'ldown'578 uu_ldown = 0.5*(rhodz(ij,l)+rhodz(ij+t_ldown,l))*u(ij+u_ldown,l) &579 -0.5*(F_el(ij+u_ldown,l)+F_el(ij+u_ldown,l+1))580 hflux(ij+u_ldown,l) = uu_ldown*le_de(ij+u_ldown)581 du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))582 END DO583 END DO584 585 #undef V_EL586 #undef G_EL587 #undef BERNI588 589 END IF ! dysl590 591 CALL trace_end("compute_caldyn_slow_NH")592 593 END SUBROUTINE compute_caldyn_slow_NH594 595 435 END MODULE caldyn_kernels_hevi_mod
Note: See TracChangeset
for help on using the changeset viewer.