source: codes/icosagcm/devel/src/unstructured/data_unstructured.F90 @ 663

MODULE data_unstructured_mod
  USE ISO_C_BINDING
  USE OMP_LIB
  IMPLICIT NONE
  SAVE

#define BINDC_(thename) BIND(C, name=#thename)
#define BINDC(thename) BINDC_(dynamico_ ## thename)

#define DBL REAL(C_DOUBLE)
#define DOUBLE1(m) DBL, DIMENSION(m)
#define DOUBLE2(m,n) DBL, DIMENSION(m,n)
#define DOUBLE3(m,n,p) DBL, DIMENSION(m,n,p)
#define DOUBLE4(m,n,p,q) DBL, DIMENSION(m,n,p,q)
#define INDEX INTEGER(C_INT)

  INTEGER, PARAMETER :: eta_mass=1, eta_lag=2, &
       thermo_theta=1, thermo_entropy=2, thermo_moist=3, thermo_boussinesq=4, &
       caldyn_vert_cons=1, max_nb_stage=5
  INTEGER(C_INT),  BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, &
       caldyn_vert_variant=caldyn_vert_cons, nb_threads=0, nb_stage=0
  LOGICAL(C_BOOL), BIND(C) :: hydrostatic=.TRUE., debug_hevi_solver=.TRUE.

  INDEX, BIND(C) :: llm, nqdyn, edge_num, primal_num, dual_num, &
       max_primal_deg, max_dual_deg, max_trisk_deg
  INDEX, ALLOCATABLE :: & ! deg(ij) = nb of vertices = nb of edges of primal/dual cell ij
       primal_deg(:), primal_edge(:,:), primal_vertex(:,:), primal_ne(:,:), & 
       dual_deg(:), dual_edge(:,:), dual_vertex(:,:), dual_ne(:,:), &
       trisk_deg(:), trisk(:,:), &
       left(:), right(:), up(:), down(:)
  ! left and right are adjacent primal cells
  ! flux is positive when going from left to right
  ! up and down are adjacent dual cells
  ! circulation is positive when going from down to up

  DBL, BIND(C) :: elapsed, g, ptop, cpp, cppv, Rd, Rv, preff, Treff, &
       pbot, Phi_bot, rho_bot
  DBL :: kappa
  DOUBLE1(max_nb_stage), BIND(C)              :: tauj       ! diagonal of fast Butcher tableau
  DOUBLE2(max_nb_stage,max_nb_stage), BIND(C) :: cslj, cflj ! slow and fast modified Butcher tableaus
  DOUBLE1(:), ALLOCATABLE            :: le_de, fv, Av, Ai
  DOUBLE2(:,:), ALLOCATABLE          :: Riv2, wee, ap,bp, mass_bl, mass_dak, mass_dbk

  INTEGER(C_INT), BIND(C) :: comm_icosa

  INTEGER, PARAMETER :: id_dev1=1, id_dev2=2, &
       id_pvort_only=3, id_slow_hydro=4, id_fast=5, id_coriolis=6, id_theta=7, id_geopot=8, id_vert=9, &
       id_solver=10, id_slow_NH=11, id_NH_geopot=12, id_vert_NH=13, id_update=14, nb_routines=14 
  DBL, PRIVATE :: start_time, time_spent(nb_routines) ! time spent in each kernel
  INTEGER, PRIVATE :: current_id, nb_calls(nb_routines), bytes(nb_routines) ! bytes read or written be each kernel
  CHARACTER(len = 10) :: id_name(nb_routines) = &
       (/'dev1      ', 'dev2      ', &
       'pvort_only', 'slow_hydro', 'fast      ', 'coriolis  ', 'theta     ', 'geopot    ', 'vert      ', &
       'solver    ', 'slow_NH   ', 'NH_geopot ', 'vert_NH   ',  'update    ' /)

CONTAINS

  !----------------------------      PROFILING      --------------------------
  
  SUBROUTINE init_trace()
    time_spent(:)=0.
    bytes(:)=0
    nb_calls(:)=0
  END SUBROUTINE init_trace

  SUBROUTINE print_trace()
    INTEGER :: id
    DBL :: total_spent
    total_spent=SUM(time_spent)
    PRINT *, '========================= Performance metrics ========================='
    PRINT *, 'Total time spent in instrumented code (seconds) :', total_spent
    PRINT *, 'Name, #calls, %time, microsec/call, MB/sec'    
    DO id=1,nb_routines
       IF(nb_calls(id)>0) PRINT *, id_name(id), nb_calls(id), INT(100.*time_spent(id)/total_spent), &
            INT(1e6*time_spent(id)/nb_calls(id)), INT(1e-6*bytes(id)/time_spent(id))
    END DO
    CALL init_trace()
  END SUBROUTINE print_trace

  SUBROUTINE enter_trace(id, nbytes)
    INTEGER :: id, nbytes
    current_id = id
    bytes(id) = bytes(id) + nbytes
    nb_calls(id)=nb_calls(id)+1
    start_time = OMP_GET_WTIME()
  END SUBROUTINE enter_trace

  SUBROUTINE exit_trace()
    DBL :: elapsed
    elapsed = OMP_GET_WTIME()-start_time
    IF(elapsed<0.) elapsed=0.
    time_spent(current_id) = time_spent(current_id) + elapsed
  END SUBROUTINE exit_trace

  !---------------------------- CONTEXT INITIALIZATION --------------------------

#define ALLOC1(v,n1) IF(ALLOCATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1))
#define ALLOC2(v,n1,n2) IF(ALLOCATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1,n2))

  SUBROUTINE init_mesh( & 
       primal_deg_, primal_edge_, primal_ne_, &
       dual_deg_, dual_edge_, dual_ne_, dual_vertex_, &
       left_, right_, up_, down_ ,&
       trisk_deg_, trisk_) BINDC(init_mesh)
    INDEX :: primal_deg_(primal_num), primal_edge_(max_primal_deg,primal_num), &
         primal_ne_(max_primal_deg,primal_num), &
         dual_deg_(dual_num), dual_edge_(max_dual_deg,dual_num), &
         dual_ne_(max_dual_deg,dual_num), &
         dual_vertex_(max_dual_deg,dual_num), &
         trisk_deg_(edge_num), trisk_(max_trisk_deg, edge_num)
    INDEX, DIMENSION(edge_num) :: left_, right_, down_, up_

    PRINT *, 'init_mesh ...'
    PRINT *, 'Primal mesh : ', primal_num, max_primal_deg
    PRINT *, 'Dual mesh   : ', dual_num, max_dual_deg
    PRINT *, 'Edge mesh   : ', edge_num, max_trisk_deg
    PRINT *, 'Vertical levels :', llm
    ALLOC1(primal_deg, primal_num)
    ALLOC2(primal_edge, max_primal_deg,primal_num)
    ALLOC2(primal_ne, max_primal_deg,primal_num)
    ALLOC1(dual_deg,dual_num)
    ALLOC2(dual_edge, max_dual_deg,dual_num)
    ALLOC2(dual_ne, max_dual_deg,dual_num)
    ALLOC2(dual_vertex, max_dual_deg,dual_num)
    ALLOC1(trisk_deg, edge_num)
    ALLOC2(trisk, max_trisk_deg, edge_num)
    PRINT *, SHAPE(trisk), edge_num
    ALLOC1(left, edge_num)
    ALLOC1(right, edge_num)
    ALLOC1(up, edge_num)
    ALLOC1(down, edge_num)
    primal_deg(:) = primal_deg_(:)
    primal_edge(:,:) = primal_edge_(:,:)
    primal_ne(:,:) = primal_ne_(:,:)
    dual_deg(:) = dual_deg_(:)
    dual_edge(:,:) = dual_edge_(:,:)
    dual_ne(:,:) = dual_ne_(:,:)
    dual_vertex(:,:) = dual_vertex_(:,:)
    IF(MINVAL(dual_deg)<3) THEN
       STOP 'At least one dual cell has less than 3 vertices'
    END IF
    IF(MINVAL(primal_deg)<3) THEN
       STOP 'At least one primal cell has less than 3 vertices'
    END IF
    left(:)=left_(:)
    right(:)=right_(:)
    down(:)=down_(:)
    up=up_(:)
    trisk_deg(:)=trisk_deg_(:)
    trisk(:,:)=trisk_(:,:)
    PRINT *, MAXVAL(primal_edge), edge_num
    PRINT *, MAXVAL(dual_edge), edge_num
    PRINT *, MAXVAL(dual_vertex), dual_num
    PRINT *, MAXVAL(trisk), edge_num
    PRINT *, MAX(MAXVAL(left),MAXVAL(right)), primal_num
    PRINT *, MAX(MAXVAL(up),MAXVAL(down)), dual_num
    PRINT *, SHAPE(trisk), edge_num
    PRINT *,' ... Done.'
  END SUBROUTINE init_mesh

  SUBROUTINE init_metric(Ai_, Av_, fv_, le_de_, Riv2_, wee_) BINDC(init_metric)
    DOUBLE1(primal_num) :: Ai_
    DOUBLE1(dual_num)   :: Av_, fv_
    DOUBLE1(edge_num)   :: le_de_
    DOUBLE2(max_dual_deg,dual_num) :: Riv2_
    DOUBLE2(max_trisk_deg,edge_num) :: wee_
    PRINT *, 'init_metric ...'
    ALLOC1(Ai,primal_num)
    ALLOC1(Av,dual_num)
    ALLOC1(fv,dual_num)
    ALLOC1(le_de,edge_num)
    ALLOC2(Riv2, max_dual_deg, dual_num)
    ALLOC2(wee, max_trisk_deg, edge_num)
    Ai(:) = Ai_(:)
    Av(:) = Av_(:)
    fv(:) = fv_(:)
    le_de(:) = le_de_(:)
    Riv2(:,:)=Riv2_(:,:)
    wee(:,:) = wee_(:,:)
    PRINT *, MAXVAL(ABS(Ai))
    PRINT *, MAXVAL(ABS(Av))
    PRINT *, MAXVAL(ABS(fv))
    PRINT *, MAXVAL(ABS(le_de))
    PRINT *, MAXVAL(ABS(Riv2))
    PRINT *, MAXVAL(ABS(wee))
    PRINT *, MINVAL(right), MAXVAL(right)
    PRINT *, MINVAL(right), MAXVAL(left)
    PRINT *,' ... Done.'
    IF(nb_threads==0) nb_threads=OMP_GET_MAX_THREADS()
    PRINT *,'OpenMP : max_threads, num_procs, nb_threads', OMP_GET_MAX_THREADS(), OMP_GET_NUM_PROCS(), nb_threads
  END SUBROUTINE init_metric
  !
  SUBROUTINE show_openmp() BINDC(show_openmp)
    PRINT *,'OpenMP : max_threads, num_procs', OMP_GET_MAX_THREADS(), OMP_GET_NUM_PROCS()
  END SUBROUTINE show_openmp
  !
  SUBROUTINE init_params() BINDC(init_params)
    PRINT *, 'Setting physical parameters ...'
    IF(hydrostatic) THEN
       PRINT *, 'Hydrostatic dynamics (HPE)'
    ELSE
       PRINT *, 'Non-hydrostatic dynamics (Euler)'
    END IF
    kappa = Rd/cpp
    PRINT *, 'g = ',g
    PRINT *, 'preff = ',preff
    PRINT *, 'Treff = ',Treff
    PRINT *, 'Rd = ',Rd
    PRINT *, 'cpp = ',cpp
    PRINT *, 'kappa = ',kappa
    PRINT *, '... Done'
  END SUBROUTINE init_params
  !
  SUBROUTINE init_hybrid(bl,dak,dbk) BINDC(init_hybrid)
    DOUBLE2(llm+1, primal_num) :: bl
    DOUBLE2(llm, primal_num) :: dak,dbk
    PRINT *, 'Setting hybrid coefficients ...'
    ALLOC2(mass_bl, llm+1, primal_num)
    ALLOC2(mass_dak, llm, primal_num)
    ALLOC2(mass_dbk, llm, primal_num)
    mass_bl(:,:)  = bl(:,:)
    mass_dak(:,:) = dak(:,:)
    mass_dbk(:,:) = dbk(:,:)
    PRINT *, '... Done, llm = ', llm
  END SUBROUTINE Init_hybrid

END MODULE data_unstructured_mod