source: codes/icosagcm/devel/src/dynamics/caldyn_hevi.f90 @ 935

MODULE caldyn_hevi_mod
  USE icosa
  USE transfert_mod
  USE caldyn_vars_mod
  USE compute_caldyn_vert_mod, ONLY : compute_caldyn_vert => compute_caldyn_vert_manual
  USE compute_caldyn_vert_NH_mod, ONLY : compute_caldyn_vert_NH => compute_caldyn_vert_NH_manual
  USE compute_caldyn_kv_mod, ONLY : compute_caldyn_kv
  USE compute_caldyn_Coriolis_mod, ONLY : compute_caldyn_Coriolis => compute_caldyn_Coriolis_manual
  USE compute_caldyn_slow_hydro_mod, ONLY : compute_caldyn_slow_hydro => compute_caldyn_slow_hydro_manual
  USE compute_caldyn_slow_NH_mod, ONLY : compute_caldyn_slow_NH
  USE compute_caldyn_solver_mod, ONLY : compute_caldyn_solver
  USE compute_NH_geopot_mod, ONLY : compute_NH_geopot
  IMPLICIT NONE
  PRIVATE
  PUBLIC caldyn_hevi

CONTAINS
  
  SUBROUTINE caldyn_hevi(write_out,tau, f_phis, f_ps, f_mass, f_theta_rhodz, f_u, f_q, &
       f_W, f_geopot, f_hflux, f_wflux, f_dps, f_dmass, f_dtheta_rhodz, &
       f_du_slow, f_du_fast, f_dPhi_slow, f_dPhi_fast, f_dW_slow, f_dW_fast) 
    USE icosa
    USE observable_mod
    USE disvert_mod, ONLY : caldyn_eta, eta_mass
    USE vorticity_mod
    USE kinetic_mod
    USE theta2theta_rhodz_mod
    USE wind_mod
    USE mpipara
    USE trace
    USE omp_para
    USE output_field_mod
    USE checksum_mod
    USE compute_caldyn_mod, ONLY : compute_pvort_only, compute_theta, &
         compute_geopot, compute_caldyn_fast
    IMPLICIT NONE
    LOGICAL,INTENT(IN)    :: write_out
    REAL(rstd), INTENT(IN) :: tau
    TYPE(t_field),POINTER :: f_phis(:)
    TYPE(t_field),POINTER :: f_ps(:)
    TYPE(t_field),POINTER :: f_mass(:)
    TYPE(t_field),POINTER :: f_theta_rhodz(:)
    TYPE(t_field),POINTER :: f_u(:)
    TYPE(t_field),POINTER :: f_q(:)
    TYPE(t_field),POINTER :: f_W(:)
    TYPE(t_field),POINTER :: f_geopot(:)
    TYPE(t_field),POINTER :: f_hflux(:), f_wflux(:)
    TYPE(t_field) :: f_dps(:)
    TYPE(t_field) :: f_dmass(:)
    TYPE(t_field) :: f_dtheta_rhodz(:)
    TYPE(t_field) :: f_du_slow(:)
    TYPE(t_field) :: f_du_fast(:)
    TYPE(t_field) :: f_dW_slow(:)
    TYPE(t_field) :: f_dW_fast(:)
    TYPE(t_field) :: f_dPhi_slow(:)
    TYPE(t_field) :: f_dPhi_fast(:)
    
    REAL(rstd),POINTER :: ps(:), dps(:), phis(:)
    REAL(rstd),POINTER :: mass(:,:), theta_rhodz(:,:,:), dtheta_rhodz(:,:,:)
    REAL(rstd),POINTER :: du(:,:), dW(:,:), dPhi(:,:), hflux(:,:), wflux(:,:)
    REAL(rstd),POINTER :: u(:,:), w(:,:), qu(:,:), qv(:,:), Kv(:,:), hv(:,:)
    REAL(rstd),POINTER :: berni(:,:), Ftheta(:,:)
! temporary shared variable
    REAL(rstd),POINTER  :: theta(:,:,:)  
    REAL(rstd),POINTER  :: pk(:,:)
    REAL(rstd),POINTER  :: geopot(:,:)
    REAL(rstd),POINTER  :: convm(:,:) 
    REAL(rstd),POINTER  :: wwuu(:,:)
    REAL(rstd),POINTER  :: F_el(:,:), gradPhi2(:,:), w_il(:,:) , W_etadot(:,:), pres(:,:), m_il(:,:)
    INTEGER :: ind
    LOGICAL,SAVE :: first=.TRUE.
!$OMP THREADPRIVATE(first)
    
    IF (first) THEN 
      first=.FALSE.
      IF(caldyn_eta==eta_mass) THEN
         CALL init_message(f_ps,req_i1,req_ps)
      ELSE
         CALL init_message(f_mass,req_i1,req_mass)
      END IF
      CALL init_message(f_theta_rhodz,req_i1,req_theta_rhodz)
      CALL init_message(f_u,req_e1_vect,req_u)
      CALL init_message(f_qu,req_e1_scal,req_qu)
      IF(caldyn_kinetic==kinetic_consistent) CALL init_message(f_Kv,req_z1_scal,req_Kv)
      IF(.NOT.hydrostatic) THEN
         CALL init_message(f_geopot,req_i1,req_geopot)
         CALL init_message(f_w,req_i1,req_w)
      END IF
    ENDIF
    
    CALL trace_start("caldyn")
    
    IF(caldyn_eta==eta_mass) THEN
       CALL send_message(f_ps,req_ps) ! COM00
       CALL wait_message(req_ps) ! COM00
    ELSE
       CALL send_message(f_mass,req_mass) ! COM00
       CALL wait_message(req_mass) ! COM00
    END IF
    CALL send_message(f_theta_rhodz,req_theta_rhodz) ! COM01
    CALL wait_message(req_theta_rhodz) ! COM01 Moved from caldyn_pvort

    IF(.NOT.hydrostatic) THEN
       CALL send_message(f_geopot,req_geopot) ! COM03
       CALL wait_message(req_geopot) ! COM03
       CALL send_message(f_w,req_w) ! COM04
       CALL wait_message(req_w) ! COM04
    END IF
    
    DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       ps=f_ps(ind)
       theta_rhodz=f_theta_rhodz(ind)
       mass=f_mass(ind)
       theta = f_theta(ind)
       CALL compute_theta(ps,theta_rhodz, mass,theta)
       pk = f_pk(ind)
       geopot = f_geopot(ind)
       du=f_du_fast(ind)
       IF(hydrostatic) THEN
          du(:,:)=0.
          CALL compute_geopot(mass,theta, ps,pk,geopot)
       ELSE
          phis = f_phis(ind)
          W = f_W(ind)
          dW = f_dW_fast(ind)
          dPhi = f_dPhi_fast(ind)
          ! reuse buffers
          m_il = f_wil(ind)
          pres = f_gradPhi2(ind)
          CALL compute_caldyn_solver(tau,phis, mass,theta,pk,geopot,W, m_il,pres, dPhi,dW,du) ! computes d(Phi,W,du)_fast and updates Phi,W
       END IF
       u=f_u(ind)
       berni=f_berni(ind) ! buffer for partial Bernoulli function
       CALL compute_caldyn_fast(tau,theta,geopot, pk,berni,du,u) ! computes du_fast and updates u
    ENDDO
    
    CALL send_message(f_u,req_u) ! COM02
    CALL wait_message(req_u)   ! COM02
    
    DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       u=f_u(ind)
       mass=f_mass(ind)
       qu=f_qu(ind)
       qv=f_qv(ind)
       hv=f_hv(ind)
!       PRINT *, 'caldyn_hevi : u, mass, qu qv, hv', &
!            SHAPE(u), SHAPE(mass), SHAPE(qu), SHAPE(qv), SHAPE(hv)
       CALL compute_pvort_only(u,mass,qu,qv,hv)
       IF(caldyn_kinetic==kinetic_consistent) THEN
          Kv=f_Kv(ind)
          CALL compute_caldyn_Kv(u,Kv)
       END IF
    ENDDO
    
    CALL send_message(f_qu,req_qu) ! COM03
    CALL wait_message(req_qu) ! COM03

    IF(caldyn_kinetic==kinetic_consistent) THEN
       CALL transfert_request(f_hv,req_z1_scal)
       CALL send_message(f_Kv,req_Kv)
       CALL wait_message(req_Kv)
    END IF

    DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       u=f_u(ind)
       mass=f_mass(ind)
       theta = f_theta(ind)
       qu=f_qu(ind)
       hflux=f_hflux(ind)
       convm = f_dmass(ind)
       dtheta_rhodz=f_dtheta_rhodz(ind)
       du=f_du_slow(ind)

       IF(hydrostatic) THEN
          hv=f_hv(ind)
          Kv=f_Kv(ind)
          berni=f_berni(ind)
          CALL compute_caldyn_slow_hydro(.TRUE., u,mass,hv,Kv, berni, hflux,du)
       ELSE
          W = f_W(ind)
          dW = f_dW_slow(ind)
          geopot = f_geopot(ind)
          dPhi = f_dPhi_slow(ind)
          F_el = f_Fel(ind)
          gradPhi2 = f_gradPhi2(ind)
          w_il = f_wil(ind)
          CALL compute_caldyn_slow_NH(u,mass,geopot,W, F_el,gradPhi2,w_il, hflux,du,dPhi,dW)
       END IF
       Ftheta=f_Ftheta(ind)
       CALL compute_caldyn_Coriolis(hflux,theta,qu, Ftheta, convm,dtheta_rhodz,du)
       IF(caldyn_eta==eta_mass) THEN
          wflux=f_wflux(ind)
          wwuu=f_wwuu(ind)
          dps=f_dps(ind)
          CALL compute_caldyn_vert(u,theta,mass,convm, wflux,wwuu, dps, dtheta_rhodz, du)
          IF(.NOT.hydrostatic) THEN
             W_etadot=f_Wetadot(ind)
             CALL compute_caldyn_vert_NH(mass,geopot,W,wflux, W_etadot, du,dPhi,dW)
          END IF
       END IF
    ENDDO
    
!$OMP BARRIER
    !    CALL check_mass_conservation(f_ps,f_dps)
    CALL trace_end("caldyn_hevi")
!!$OMP BARRIER
    
  END SUBROUTINE caldyn_hevi

END MODULE caldyn_hevi_mod