Changeset 734 for codes


Ignore:
Timestamp:
08/27/18 13:51:26 (6 years ago)
Author:
dubos
Message:

devel : consistent discretization of kinetic energy

Location:
codes/icosagcm/devel/src/dynamics
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90

    r733 r734  
    5151    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:) 
    5252    REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:) 
    53     REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:) 
     53    REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:), Kv(:,:) 
    5454 
    5555! temporary shared variable 
     
    7474      CALL init_message(f_u,req_e1_vect,req_u) 
    7575      CALL init_message(f_qu,req_e1_scal,req_qu) 
     76      IF(caldyn_kinetic==kinetic_consistent) CALL init_message(f_Kv,req_z1_scal,req_Kv) 
    7677      IF(.NOT.hydrostatic) THEN 
    7778         CALL init_message(f_geopot,req_i1,req_geopot) 
     
    140141       qv=f_qv(ind) 
    141142       CALL compute_pvort_only(u,mass,qu,qv) 
     143       IF(caldyn_kinetic==kinetic_consistent) THEN 
     144          Kv=f_Kv(ind) 
     145          CALL compute_caldyn_Kv(u,Kv) 
     146       END IF 
    142147    ENDDO 
    143148     
    144149    CALL send_message(f_qu,req_qu) ! COM03 
    145150    CALL wait_message(req_qu) ! COM03 
    146         
     151 
     152    IF(caldyn_kinetic==kinetic_consistent) THEN 
     153       CALL send_message(f_Kv,req_Kv) 
     154       CALL wait_message(req_Kv) 
     155    END IF 
     156 
    147157    DO ind=1,ndomain 
    148158       IF (.NOT. assigned_domain(ind)) CYCLE 
     
    159169 
    160170       IF(hydrostatic) THEN 
    161           CALL compute_caldyn_slow_hydro(u,mass,hflux,du, .TRUE.) 
     171          Kv=f_Kv(ind) 
     172          CALL compute_caldyn_slow_hydro(u,mass,hflux,Kv,du, .TRUE.) 
    162173       ELSE 
    163174          W = f_W(ind) 
  • codes/icosagcm/devel/src/dynamics/caldyn_kernels_hevi.F90

    r733 r734  
    1313  LOGICAL, SAVE :: debug_hevi_solver = .FALSE. 
    1414 
    15   PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Coriolis, & 
     15  PUBLIC :: compute_theta, compute_pvort_only, compute_caldyn_Kv, compute_caldyn_Coriolis, & 
    1616       compute_caldyn_slow_hydro, compute_caldyn_slow_NH, & 
    1717       compute_caldyn_solver, compute_caldyn_fast 
     
    9999 
    100100  END SUBROUTINE compute_pvort_only 
     101 
     102  SUBROUTINE compute_caldyn_kv(ue, Kv) 
     103    REAL(rstd),INTENT(IN) :: ue(3*iim*jjm,llm) 
     104    REAL(rstd),INTENT(OUT) :: Kv(2*iim*jjm,llm) 
     105    REAL(rstd) :: ue2(3*iim*jjm), dem2(3*iim*jjm), r2_Av(2*iim*jjm), rad2 
     106    INTEGER :: ij,l, u_up, u_down 
     107     
     108    u_up    = t_lup + u_right 
     109    u_down  = t_rdown + u_left 
     110     
     111    rad2=radius**2 
     112 
     113    !DIR$ SIMD 
     114    DO ij=ij_begin_ext,ij_end_ext 
     115       dem2(ij+u_right) = de(ij+u_right)**(-2) 
     116       dem2(ij+u_lup)   = de(ij+u_lup)**(-2) 
     117       dem2(ij+u_ldown) = de(ij+u_ldown)**(-2) 
     118       r2_Av(ij+z_up)   = rad2*(1./Av(ij+z_up)) 
     119       r2_Av(ij+z_down) = rad2*(1./Av(ij+z_down)) 
     120    END DO 
     121 
     122    DO l=ll_begin,ll_end 
     123       ! compute squared normal component from 1-form 
     124       !DIR$ SIMD 
     125       DO ij=ij_begin_ext,ij_end_ext 
     126          ue2(ij+u_right) = dem2(ij+u_right)* (ue(ij+u_right,l)**2) 
     127          ue2(ij+u_lup)   = dem2(ij+u_lup)  * (ue(ij+u_lup,l)**2) 
     128          ue2(ij+u_ldown) = dem2(ij+u_ldown)* (ue(ij+u_ldown,l)**2) 
     129       END DO 
     130       ! average squared normal component to vertices 
     131       !DIR$ SIMD 
     132       DO ij=ij_begin_ext,ij_end_ext 
     133          Kv(ij+z_up,l) = r2_Av(ij+z_up)*(                       & 
     134                               S1(ij,vup)*ue2(ij+u_rup) +        & 
     135                               S2(ij,vup)*ue2(ij+u_lup) +        & 
     136                               S2(ij+t_lup,vrdown)*ue2(ij+u_up)) 
     137 
     138          Kv(ij+z_down,l) = r2_Av(ij+z_down)*(                   & 
     139                               S1(ij,vdown)*ue2(ij+u_ldown) +    & 
     140                               S2(ij,vdown)*ue2(ij+u_rdown) +    & 
     141                               S2(ij+t_rdown,vlup)*ue2(ij+u_down) ) 
     142       ENDDO 
     143    ENDDO 
     144  END SUBROUTINE compute_caldyn_kv 
    101145 
    102146  SUBROUTINE compute_NH_geopot(tau, phis, m_ik, m_il, theta, W_il, Phi_il) 
     
    663707                  wee(ij+u_ldown,5,2)*hflux(ij+t_ldown+u_right,l) 
    664708 
    665              du(ij+u_right,l) = du(ij+u_right,l) + .5*uu_right 
    666              du(ij+u_lup,l)   = du(ij+u_lup,l)   + .5*uu_lup 
    667              du(ij+u_ldown,l) = du(ij+u_ldown,l)   + .5*uu_ldown 
     709             du(ij+u_right,l) = du(ij+u_right,l) + uu_right*qu(ij+u_right,l) 
     710             du(ij+u_lup,l)   = du(ij+u_lup,l)   + uu_lup*qu(ij+u_lup,l)      
     711             du(ij+u_ldown,l) = du(ij+u_ldown,l) + uu_ldown*qu(ij+u_ldown,l)  
    668712          END DO 
    669713       END DO 
     
    679723  END SUBROUTINE compute_caldyn_Coriolis 
    680724 
    681   SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,du, zero) 
     725  SUBROUTINE compute_caldyn_slow_hydro(u,rhodz,hflux,Kv,du, zero) 
    682726    LOGICAL, INTENT(IN) :: zero 
    683727    REAL(rstd),INTENT(IN)  :: u(3*iim*jjm,llm)    ! prognostic "velocity" 
     728    REAL(rstd),INTENT(IN)  :: Kv(2*iim*jjm,llm)   ! kinetic energy at vertices 
    684729    REAL(rstd),INTENT(IN)  :: rhodz(iim*jjm,llm) 
    685730    REAL(rstd),INTENT(OUT) :: hflux(3*iim*jjm,llm) ! hflux in kg/s 
     
    719764       ENDDO 
    720765       ! Compute Bernoulli=kinetic energy 
    721        !DIR$ SIMD 
    722        DO ij=ij_begin,ij_end          
    723           BERNI(ij) = & 
    724                1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    & 
    725                le_de(ij+u_rup)*u(ij+u_rup,l)**2 +          & 
    726                le_de(ij+u_lup)*u(ij+u_lup,l)**2 +          & 
    727                le_de(ij+u_left)*u(ij+u_left,l)**2 +       & 
    728                le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
    729                le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
    730        ENDDO 
     766       IF(caldyn_kinetic==kinetic_trisk) THEN 
     767          !DIR$ SIMD 
     768          DO ij=ij_begin,ij_end          
     769             BERNI(ij) = & 
     770                  1/(4*Ai(ij))*(le_de(ij+u_right)*u(ij+u_right,l)**2 +    & 
     771                                le_de(ij+u_rup)*u(ij+u_rup,l)**2     +    & 
     772                                le_de(ij+u_lup)*u(ij+u_lup,l)**2     +    & 
     773                                le_de(ij+u_left)*u(ij+u_left,l)**2   +    & 
     774                                le_de(ij+u_ldown)*u(ij+u_ldown,l)**2 +    & 
     775                                le_de(ij+u_rdown)*u(ij+u_rdown,l)**2 )   
     776          ENDDO 
     777       ELSE 
     778          !DIR$ SIMD 
     779          DO ij=ij_begin,ij_end 
     780             BERNI(ij) = Riv(ij,vup)   *Kv(ij+z_up,l)    + & 
     781                         Riv(ij,vlup)  *Kv(ij+z_lup,l)   + & 
     782                         Riv(ij,vldown)*Kv(ij+z_ldown,l) + & 
     783                         Riv(ij,vdown) *Kv(ij+z_down,l)  + & 
     784                         Riv(ij,vrdown)*Kv(ij+z_rdown,l) + & 
     785                         Riv(ij,vrup)  *Kv(ij+z_rup,l) 
     786          END DO 
     787       END IF 
    731788       ! Compute du=-grad(Bernoulli) 
    732789       IF(zero) THEN 
  • codes/icosagcm/devel/src/dynamics/caldyn_vars.f90

    r733 r734  
    2626                           f_Fel(:), f_gradPhi2(:), f_wil(:), f_Wetadot(:) 
    2727 
    28   TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu, req_geopot, req_w 
     28  TYPE(t_message) :: req_ps, req_mass, req_theta_rhodz, req_u, req_qu, req_geopot, req_w, req_Kv 
    2929 
    3030END MODULE caldyn_vars_mod 
Note: See TracChangeset for help on using the changeset viewer.