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

Last change on this file since 836 was 836, checked in by dubos, 5 years ago

devel : Cp(T) thermodynamics (TBC)
File size: 10.7 KB
Line 
1MODULE data_unstructured_mod
2  USE ISO_C_BINDING
3  USE earth_const, ONLY : thermo_theta
4  USE mpipara, ONLY : is_mpi_master
5  USE grid_param, ONLY : llm, nqdyn
6#ifdef CPP_USING_OMP
7  USE OMP_LIB
8#endif
9  IMPLICIT NONE
10  SAVE
11
12
13#include "unstructured.h90"
14
15  INTEGER, PARAMETER :: eta_mass=1, eta_lag=2, &
16!       thermo_theta=1, thermo_entropy=2, thermo_moist=3, thermo_boussinesq=4, &
17       caldyn_vert_cons=1, max_nb_stage=5
18  INDEX,  BIND(C) :: caldyn_eta=eta_lag, &
19       caldyn_vert_variant=caldyn_vert_cons, nb_threads=0, nb_stage=0
20!  INDEX,  BIND(C) :: caldyn_thermo=thermo_theta, caldyn_eta=eta_lag, &
21!       caldyn_vert_variant=caldyn_vert_cons, nb_threads=0, nb_stage=0
22  LOGICAL(C_BOOL), BIND(C) :: hydrostatic=.TRUE., debug_on=.FALSE.
23  LOGICAL(C_BOOL), BIND(C, NAME='debug_hevi_solver') :: debug_hevi_solver_=.TRUE.
24
25#ifdef CPP_MIXED_PREC
26  LOGICAL(C_BOOL), BIND(C) :: mixed_precision=.TRUE.
27#else
28  LOGICAL(C_BOOL), BIND(C) :: mixed_precision=.FALSE.
29#endif
30
31  INDEX, BIND(C) :: edge_num, primal_num, dual_num, &
32       max_primal_deg, max_dual_deg, max_trisk_deg
33  INDEX, ALLOCATABLE :: & ! deg(ij) = nb of vertices = nb of edges of primal/dual cell ij
34       primal_deg(:), primal_edge(:,:), primal_vertex(:,:), primal_ne(:,:), & 
35       dual_deg(:), dual_edge(:,:), dual_vertex(:,:), dual_ne(:,:), &
36       trisk_deg(:), trisk(:,:), &
37       left(:), right(:), up(:), down(:)
38  ! left and right are adjacent primal cells
39  ! flux is positive when going from left to right
40  ! up and down are adjacent dual cells
41  ! circulation is positive when going from down to up
42
43  TIME, PARAMETER :: print_trace_interval = 1.
44  TIME, BIND(C) :: elapsed
45  NUM, BIND(C) :: ptop, pbot, Phi_bot, rho_bot
46!  NUM, BIND(C) :: g, ptop, cpp, cppv, Rd, Rv, preff, Treff, pbot, Phi_bot, rho_bot
47!  NUM :: kappa
48  NUM1(max_nb_stage), BIND(C)              :: tauj       ! diagonal of fast Butcher tableau
49  NUM2(max_nb_stage,max_nb_stage), BIND(C) :: cslj, cflj ! slow and fast modified Butcher tableaus
50  NUM1(:), ALLOCATABLE            :: le, le_de, fv, Av, Ai
51  NUM2(:,:), ALLOCATABLE          :: centroid, xyz_v, Riv2, wee, ap,bp, mass_bl, mass_dak, mass_dbk
52
53  INTEGER(C_INT), BIND(C) :: comm_icosa
54
55  INTEGER, PARAMETER :: id_dev1=1, id_dev2=2, &
56       id_pvort_only=3, id_slow_hydro=4, id_fast=5, id_coriolis=6, id_theta=7, id_geopot=8, id_vert=9, &
57       id_solver=10, id_slow_NH=11, id_NH_geopot=12, id_vert_NH=13, id_update=14, id_halo=15, &
58       id_scalar_laplacian=16, nb_routines=16 
59  TIME, PRIVATE :: start_time, time_spent(nb_routines) ! time spent in each kernel
60  INTEGER, PRIVATE :: current_id, nb_calls(nb_routines)
61  INTEGER(KIND=8), PRIVATE :: bytes(nb_routines) ! bytes read or written by each kernel
62  CHARACTER(len = 10) :: id_name(nb_routines) = &
63       (/'dev1      ', 'dev2      ', &
64       'pvort_only', 'slow_hydro', 'fast      ', 'coriolis  ', 'theta     ', 'geopot    ', 'vert      ', &
65       'solver    ', 'slow_NH   ', 'NH_geopot ', 'vert_NH   ',  'update    ', 'halo_xchg ', 'scalar_lap' /)
66
67  INTEGER, PARAMETER ::transfer_primal=1, transfer_edge=2, transfer_dual=3, transfer_max=3
68  TYPE Halo_transfer
69     INTEGER :: ranks ! size of arrays rank, len
70     INTEGER, ALLOCATABLE :: rank(:), & ! MPI ranks to communicate with
71          num(:), & ! number of cells to send to / receive from other MPI ranks
72          cells(:) ! local indices of cells to send/receive
73     NUM, ALLOCATABLE :: buf2(:,:)
74  END TYPE Halo_transfer
75  TYPE(Halo_transfer), TARGET :: send_info(transfer_max), recv_info(transfer_max)
76
77CONTAINS
78
79  !----------------------------      PROFILING      --------------------------
80
81#ifndef CPP_USING_OMP
82  FUNCTION omp_get_wtime()
83    TIME :: omp_get_wtime
84    CALL CPU_TIME(omp_get_wtime)
85  END FUNCTION omp_get_wtime
86
87  FUNCTION omp_get_num_procs()
88    INTEGER :: omp_get_num_procs
89    omp_get_num_procs=1
90  END FUNCTION omp_get_num_procs
91
92  FUNCTION omp_get_max_threads()
93    INTEGER :: omp_get_max_threads
94    omp_get_max_threads=1
95  END FUNCTION omp_get_max_threads
96#endif
97
98  SUBROUTINE init_trace()
99    !$OMP MASTER
100    time_spent(:)=0.
101    bytes(:)=0
102    nb_calls(:)=0
103    !$OMP END MASTER
104  END SUBROUTINE init_trace
105
106  SUBROUTINE print_trace_() BIND(C, name='dynamico_print_trace')
107    INTEGER :: id
108    TIME :: total_spent
109    total_spent=SUM(time_spent)
110    IF(is_mpi_master) THEN
111       PRINT *, '========================= Performance metrics ========================='
112       PRINT *, 'Total time spent in instrumented code (seconds) :', total_spent
113       PRINT *, 'Name, #calls, %time, microsec/call, MB/sec'   
114       DO id=1,nb_routines
115          IF(nb_calls(id)>0) PRINT *, id_name(id), nb_calls(id), INT(100.*time_spent(id)/total_spent), &
116               INT(1e6*time_spent(id)/nb_calls(id)), INT(1e-6*bytes(id)/time_spent(id))
117       END DO
118    END IF
119  END SUBROUTINE print_trace_
120
121  SUBROUTINE print_trace()
122    !$OMP MASTER
123       IF(SUM(time_spent)>print_trace_interval) THEN
124          CALL print_trace_
125          CALL init_trace()
126       END IF
127    !$OMP END MASTER
128  END SUBROUTINE print_trace
129
130  SUBROUTINE enter_trace(id, nbytes)
131    INTEGER :: id, nbytes
132    !$OMP MASTER
133    current_id = id
134    bytes(id) = bytes(id) + nbytes
135    nb_calls(id)=nb_calls(id)+1
136    start_time = OMP_GET_WTIME()
137    !$OMP END MASTER
138  END SUBROUTINE enter_trace
139
140  SUBROUTINE exit_trace()
141    TIME :: elapsed
142    !$OMP MASTER
143    elapsed = OMP_GET_WTIME()-start_time
144    IF(elapsed<0.) elapsed=0.
145    time_spent(current_id) = time_spent(current_id) + elapsed
146    !$OMP END MASTER
147  END SUBROUTINE exit_trace
148
149  !---------------------------- CONTEXT INITIALIZATION --------------------------
150
151#define ALLOC1(v,n1) IF(ALLOCATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1))
152#define ALLOC2(v,n1,n2) IF(ALLOCATED(v)) DEALLOCATE(v) ; ALLOCATE(v(n1,n2))
153
154  SUBROUTINE init_mesh( & 
155       primal_deg_, primal_edge_, primal_ne_, &
156       dual_deg_, dual_edge_, dual_ne_, dual_vertex_, &
157       left_, right_, up_, down_ ,&
158       trisk_deg_, trisk_) BINDC(init_mesh)
159    INDEX :: primal_deg_(primal_num), primal_edge_(max_primal_deg,primal_num), &
160         primal_ne_(max_primal_deg,primal_num), &
161         dual_deg_(dual_num), dual_edge_(max_dual_deg,dual_num), &
162         dual_ne_(max_dual_deg,dual_num), &
163         dual_vertex_(max_dual_deg,dual_num), &
164         trisk_deg_(edge_num), trisk_(max_trisk_deg, edge_num)
165    INDEX, DIMENSION(edge_num) :: left_, right_, down_, up_
166
167    IF(is_mpi_master) THEN
168       PRINT *, 'init_mesh ...'
169       PRINT *, 'Primal mesh : ', primal_num, max_primal_deg
170       PRINT *, 'Dual mesh   : ', dual_num, max_dual_deg
171       PRINT *, '       Edge mesh   : ', edge_num, max_trisk_deg
172       PRINT *, 'Vertical levels :', llm
173    END IF
174    ALLOC1(primal_deg, primal_num)
175    ALLOC2(primal_edge, max_primal_deg,primal_num)
176    ALLOC2(primal_ne, max_primal_deg,primal_num)
177    ALLOC1(dual_deg,dual_num)
178    ALLOC2(dual_edge, max_dual_deg,dual_num)
179    ALLOC2(dual_ne, max_dual_deg,dual_num)
180    ALLOC2(dual_vertex, max_dual_deg,dual_num)
181    ALLOC1(trisk_deg, edge_num)
182    ALLOC2(trisk, max_trisk_deg, edge_num)
183    ALLOC1(left, edge_num)
184    ALLOC1(right, edge_num)
185    ALLOC1(up, edge_num)
186    ALLOC1(down, edge_num)
187    primal_deg(:) = primal_deg_(:)
188    primal_edge(:,:) = primal_edge_(:,:)
189    primal_ne(:,:) = primal_ne_(:,:)
190    dual_deg(:) = dual_deg_(:)
191    dual_edge(:,:) = dual_edge_(:,:)
192    dual_ne(:,:) = dual_ne_(:,:)
193    dual_vertex(:,:) = dual_vertex_(:,:)
194    IF(MINVAL(dual_deg)<2) THEN
195       STOP 'At least one dual cell has less than 2 vertices'
196    END IF
197    IF(MINVAL(primal_deg)<2) THEN
198       STOP 'At least one primal cell has less than 2 vertices'
199    END IF
200    left(:)=left_(:)
201    right(:)=right_(:)
202    down(:)=down_(:)
203    up=up_(:)
204    trisk_deg(:)=trisk_deg_(:)
205    trisk(:,:)=trisk_(:,:)
206    IF(is_mpi_master) THEN
207       PRINT *, MAXVAL(primal_edge), edge_num
208       PRINT *, MAXVAL(dual_edge), edge_num
209       PRINT *, MAXVAL(dual_vertex), dual_num
210       PRINT *, MAXVAL(trisk), edge_num
211       PRINT *, MAX(MAXVAL(left),MAXVAL(right)), primal_num
212       PRINT *, MAX(MAXVAL(up),MAXVAL(down)), dual_num
213       PRINT *, SHAPE(trisk), edge_num
214       PRINT *,' ... Done.'
215    END IF
216  END SUBROUTINE init_mesh
217
218  ! Input arrays to init_metric and init_hybrid are declared DBL
219  ! => always float64 on the Python side
220  ! They are copied to Fortran arrays of type NUM (float or double)
221
222  SUBROUTINE init_metric(Ai_, Av_, fv_, le_de_, Riv2_, wee_) BINDC(init_metric)
223    DBL :: Ai_(primal_num), Av_(dual_num), fv_(dual_num), le_de_(edge_num), &
224         Riv2_(max_dual_deg,dual_num), wee_(max_trisk_deg,edge_num)
225    IF(is_mpi_master) PRINT *, 'init_metric ...'
226    ALLOC1(Ai,primal_num)
227    ALLOC1(Av,dual_num)
228    ALLOC1(fv,dual_num)
229    ALLOC1(le_de,edge_num)
230    ALLOC2(Riv2, max_dual_deg, dual_num)
231    ALLOC2(wee, max_trisk_deg, edge_num)
232    Ai(:) = Ai_(:)
233    Av(:) = Av_(:)
234    fv(:) = fv_(:)
235    le_de(:) = le_de_(:)
236    Riv2(:,:)=Riv2_(:,:)
237    wee(:,:) = wee_(:,:)
238    IF(is_mpi_master) THEN
239       PRINT *, 'Max Ai : ',    MAXVAL(ABS(Ai))
240       PRINT *, 'Max Av : ',    MAXVAL(ABS(Av))
241       PRINT *, 'Max fv : ',    MAXVAL(ABS(fv))
242       PRINT *, 'Max le_de : ', MAXVAL(ABS(le_de))
243       PRINT *, 'Max Riv2 : ',  MAXVAL(ABS(Riv2))
244       PRINT *, 'Max wee : ',   MAXVAL(ABS(wee))
245       PRINT *, MINVAL(right),  MAXVAL(right)
246       PRINT *, MINVAL(right),  MAXVAL(left)
247       PRINT *,' ... Done.'
248       IF(nb_threads==0) nb_threads=OMP_GET_MAX_THREADS()
249       PRINT *,'OpenMP : max_threads, num_procs, nb_threads', OMP_GET_MAX_THREADS(), OMP_GET_NUM_PROCS(), nb_threads
250    END IF
251  END SUBROUTINE init_metric
252  !
253  SUBROUTINE show_openmp() BINDC(show_openmp)
254    PRINT *,'OpenMP : max_threads, num_procs', OMP_GET_MAX_THREADS(), OMP_GET_NUM_PROCS()
255  END SUBROUTINE show_openmp
256  !
257  SUBROUTINE init_params() BINDC(init_params)
258    USE earth_const
259    kappa = Rd/cpp
260    IF(is_mpi_master) THEN
261       PRINT *, 'Setting physical parameters ...'
262       IF(hydrostatic) THEN
263          PRINT *, 'Hydrostatic dynamics (HPE)'
264       ELSE
265          PRINT *, 'Non-hydrostatic dynamics (Euler)'
266       END IF
267       PRINT *, 'g = ',g
268       PRINT *, 'preff = ',preff
269       PRINT *, 'Treff = ',Treff
270       PRINT *, 'Rd = ',Rd
271       PRINT *, 'cpp = ',cpp
272       PRINT *, 'kappa = ',kappa
273       PRINT *, '... Done'
274    END IF
275    CALL init_trace
276  END SUBROUTINE init_params
277  !
278  SUBROUTINE init_hybrid(bl,dak,dbk) BINDC(init_hybrid)
279    DBL :: bl(llm+1, primal_num), &
280         dak(llm, primal_num), dbk(llm, primal_num)
281    IF(is_mpi_master) PRINT *, 'Setting hybrid coefficients ...'
282    ALLOC2(mass_bl, llm+1, primal_num)
283    ALLOC2(mass_dak, llm, primal_num)
284    ALLOC2(mass_dbk, llm, primal_num)
285    mass_bl(:,:)  = bl(:,:)
286    mass_dak(:,:) = dak(:,:)
287    mass_dbk(:,:) = dbk(:,:)
288    IF(is_mpi_master) PRINT *, '... Done, llm = ', llm
289  END SUBROUTINE Init_hybrid
290
291END MODULE data_unstructured_mod
Note: See TracBrowser for help on using the repository browser.