Context Navigation

source: codes/icosagcm/devel/src/unstructured/timestep_unstructured.F90 @ 683

Last change on this file since 683 was 683, checked in by dubos, 6 years ago
devel/unstructured : cleanup
File size: 12.1 KB

Line
1	MODULE timestep_unstructured_mod
2	USE ISO_C_BINDING
3	USE OMP_LIB
4	#ifdef CPP_USING_XIOS
5	USE xios
6	#endif
7	USE caldyn_unstructured_mod
8	USE transfer_unstructured_mod
9	IMPLICIT NONE
10	PRIVATE
11	SAVE
12
13	#define BINDC_(thename) BIND(C, name=#thename)
14	#define BINDC(thename) BINDC_(dynamico_ ## thename)
15
16	#define DBL REAL(C_DOUBLE)
17	#define DOUBLE1(m) DBL, DIMENSION(m)
18	#define DOUBLE2(m,n) DBL, DIMENSION(m,n)
19	#define DOUBLE3(m,n,p) DBL, DIMENSION(m,n,p)
20	#define DOUBLE4(m,n,p,q) DBL, DIMENSION(m,n,p,q)
21	#define INDEX INTEGER(C_INT)
22
23	#ifdef CPP_USING_XIOS
24	TYPE(xios_context) :: ctx_hdl
25	#endif
26
27	CONTAINS
28
29	#define FIELD_PS DOUBLE1(primal_num)
30	#define FIELD_MASS DOUBLE2(llm, primal_num)
31	#define FIELD_Z DOUBLE2(llm, dual_num)
32	#define FIELD_U DOUBLE2(llm, edge_num)
33	#define FIELD_UL DOUBLE2(llm+1, edge_num)
34	#define FIELD_THETA DOUBLE3(llm, primal_num, nqdyn)
35	#define FIELD_GEOPOT DOUBLE2(llm+1, primal_num)
36
37	#define HASNAN(field) (ANY(.NOT.ABS(field)<1e20))
38
39	SUBROUTINE caldyn_unstructured(tau, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! IN : flow state
40	theta,ps,pk,hflux,qv, & ! OUT : diags (except surface geopot : IN)
41	dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, &
42	dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(caldyn_unstructured) ! OUT : tendencies
43	DBL, VALUE :: tau
44	FIELD_MASS :: rhodz, drhodz, pk, berni ! IN, OUT, DIAG
45	FIELD_THETA :: theta_rhodz, dtheta_rhodz, theta ! IN, OUT, DIAG
46	FIELD_GEOPOT :: wflux, w, geopot, & ! DIAG, INOUT
47	dPhi_fast, dPhi_slow, dW_fast, dW_slow ! OUT
48	FIELD_U :: u,du_fast,du_slow,hflux,qu ! INOUT,OUT,OUT,DIAG
49	FIELD_Z :: qv ! DIAG
50	FIELD_PS :: ps,dmass_col,mass_col ! OUT,OUT,IN (if eta_mass) or OUT,UNUSED,UNUSED (if eta_lag)
51	! buffers for fields that need to be shared among OpenMP threads
52	! vertical size is allocated as llm+1 even if only llm is needed
53	FIELD_GEOPOT :: bufm1, bufm2, bufm3
54	FIELD_UL :: bufu1, bufu2, bufu3,bufu4
55
56	DBL :: time1,time2
57	INTEGER :: ij
58
59	! CALL CPU_TIME(time1)
60	time1=OMP_GET_WTIME()
61
62	IF(hydrostatic) THEN
63
64	!$OMP PARALLEL NUM_THREADS(nb_threads)
65	!$OMP DO SCHEDULE(STATIC)
66	DO ij=1,edge_num
67	du_fast(:,ij)=0.
68	du_slow(:,ij)=0.
69	END DO
70	!$OMP END DO
71	CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
72	CALL compute_geopot(rhodz,theta, ps,pk,geopot)
73
74	CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
75
76	CALL compute_pvort_only(rhodz,u,qv,qu)
77	CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow)
78	CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow)
79	IF(caldyn_eta == eta_mass) THEN
80	CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1)
81	END IF
82	!$OMP END PARALLEL
83
84	ELSE ! NH
85	DO ij=1,edge_num
86	du_fast(:,ij)=0.
87	du_slow(:,ij)=0.
88	END DO
89	DO ij=1,primal_num
90	wflux(1,ij)=0.
91	wflux(llm+1,ij)=0.
92	END DO
93	CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
94	CALL compute_caldyn_solver(tau,rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W,dPhi_fast,dW_fast,du_fast)
95	CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
96	CALL compute_pvort_only(rhodz,u,qv,qu)
97	CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3,bufu1,bufu2,bufu3,bufu4, hflux,du_slow,dPhi_slow,dW_slow)
98	CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow)
99	IF(caldyn_eta == eta_mass) THEN
100	CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1)
101	CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, du_slow,dPhi_slow,dW_slow)
102	END IF
103	END IF
104
105	time2=OMP_GET_WTIME()
106	! CALL CPU_TIME(time2)
107	IF(time2>time1) elapsed = elapsed + time2-time1
108
109	CALL print_trace()
110	END SUBROUTINE caldyn_unstructured
111	!
112	!----------------------------- Time stepping -------------------------------
113
114	!
115	SUBROUTINE ARK_step(nstep, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! INOUT : flow state
116	theta,ps,pk,hflux,qv, & ! OUT : diags (except surface geopot : IN)
117	dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, &
118	dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(ARK_step) ! OUT : tendencies
119	INTEGER, VALUE :: nstep ! advance by nstep time steps
120	FIELD_PS :: mass_col, ps ! OUT,IN (if eta_mass) or OUT,UNUSED (if eta_lag)
121	FIELD_MASS :: rhodz, pk, berni ! IN, DIAG
122	FIELD_THETA :: theta_rhodz, theta ! IN, DIAG
123	FIELD_U :: u,hflux,qu ! INOUT,DIAG*2
124	FIELD_GEOPOT :: geopot, w, wflux ! IN, INOUT, DIAG
125	FIELD_Z :: qv ! DIAG
126	DOUBLE2( primal_num, max_nb_stage) :: dmass_col ! OUT
127	DOUBLE3(llm, primal_num, max_nb_stage) :: drhodz ! OUT
128	DOUBLE4(llm, primal_num, nqdyn, max_nb_stage) :: dtheta_rhodz ! OUT
129	DOUBLE3(llm, edge_num, max_nb_stage) :: du_fast,du_slow ! OUT
130	DOUBLE3(llm+1, primal_num, max_nb_stage) :: &
131	dPhi_fast, dPhi_slow, dW_fast, dW_slow ! OUT
132	! buffers for fields that need to be shared among OpenMP threads
133	! vertical size is allocated as llm+1 even if only llm is needed
134	FIELD_GEOPOT :: bufm1, bufm2, bufm3
135	FIELD_UL :: bufu1, bufu2, bufu3,bufu4
136	DBL :: time1,time2
137	INTEGER :: step, stage, iq, ij
138
139	!CALL CPU_TIME(time1)
140	time1=OMP_GET_WTIME()
141
142	!$OMP PARALLEL NUM_THREADS(nb_threads)
143
144	!$OMP DO SCHEDULE(STATIC)
145	DO ij=1,primal_num
146	wflux(1,ij)=0.
147	wflux(llm+1,ij)=0.
148	END DO
149	!$OMP END DO
150
151	DO step=1, nstep
152	DO stage=1, nb_stage
153
154	CALL update_halo(transfer_primal, rhodz)
155	DO iq=1, nqdyn
156	CALL update_halo(transfer_primal, theta_rhodz(:,:,iq))
157	END DO
158
159	IF(hydrostatic) THEN
160
161	!$OMP DO SCHEDULE(STATIC)
162	DO ij=1,edge_num
163	du_fast(:,ij,stage)=0.
164	du_slow(:,ij,stage)=0.
165	END DO
166
167	CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
168	CALL compute_geopot(rhodz,theta, ps,pk,geopot)
169
170	CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u)
171	CALL update_halo(transfer_edge, u)
172	CALL compute_pvort_only(rhodz,u,qv,qu)
173	CALL update_halo(transfer_edge, qu)
174
175	CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage))
176	CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz(:,:,stage), &
177	dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage))
178	IF(caldyn_eta == eta_mass) THEN
179	CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
180	dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),bufu1)
181	END IF
182
183	ELSE ! NH
184	CALL update_halo(transfer_primal, geopot)
185	CALL update_halo(transfer_primal, W)
186
187	CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
188	CALL compute_caldyn_solver(tauj(stage),rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W, &
189	dPhi_fast(:,:,stage), dW_fast(:,:,stage), du_fast(:,:,stage))
190
191	CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage),u)
192	CALL update_halo(transfer_edge, u)
193	CALL compute_pvort_only(rhodz,u,qv,qu)
194	CALL update_halo(transfer_edge, qu)
195
196	CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3, &
197	bufu1,bufu2,bufu3,bufu4, hflux, &
198	du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage))
199	CALL compute_coriolis(hflux,theta,qu, bufu1, &
200	drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage))
201	IF(caldyn_eta == eta_mass) THEN
202	CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
203	dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), bufu1)
204	CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, &
205	du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) )
206	END IF
207	END IF ! NH
208
209	! FIXME : mass_col is computed from rhodz
210	! so that the DOFs are the same whatever caldyn_eta
211	! in DYNAMICO mass_col is prognosed rather than rhodz
212
213	#define UPDATE(clj,field,dfield) update(stage,SIZE(field),clj,field,dfield)
214
215	CALL UPDATE(cslj, rhodz, drhodz)
216	CALL UPDATE(cslj, theta_rhodz, dtheta_rhodz)
217	CALL UPDATE(cslj, u, du_slow)
218	CALL UPDATE(cflj, u, du_fast)
219
220	IF(.NOT.hydrostatic) THEN
221	CALL UPDATE(cslj, geopot, dPhi_slow)
222	CALL UPDATE(cflj, geopot, dPhi_fast)
223	CALL UPDATE(cslj, W, dW_slow)
224	CALL UPDATE(cflj, W, dW_fast)
225	END IF
226
227	END DO
228	END DO
229	!$OMP END PARALLEL
230
231	time2=OMP_GET_WTIME()
232	! CALL CPU_TIME(time2)
233	IF(time2>time1) elapsed = elapsed + time2-time1
234
235	CALL print_trace()
236
237	END SUBROUTINE ARK_step
238
239
240	SUBROUTINE update(j,sz,clj,field,dfield)
241	INTEGER :: j, sz ! stage in ARK scheme, field size
242	DOUBLE2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau
243	DOUBLE1(sz) :: field
244	DOUBLE2(sz, max_nb_stage) :: dfield
245	!
246	INTEGER :: l, ij
247	CALL enter_trace(id_update, 8sz(j+1) )
248	! PRINT *, clj(1:j,j)
249	SELECT CASE(j)
250	CASE(1)
251	!$OMP DO SCHEDULE(static)
252	DO ij=1,sz
253	field(ij) = field(ij) &
254	+ clj(1,j)*dfield(ij,1)
255	END DO
256	CASE(2)
257	!$OMP DO SCHEDULE(static)
258	DO ij=1,sz
259	field(ij) = field(ij) &
260	+ clj(1,j)*dfield(ij,1) &
261	+ clj(2,j)*dfield(ij,2)
262	END DO
263	CASE(3)
264	!$OMP DO SCHEDULE(static)
265	DO ij=1,sz
266	field(ij) = field(ij) &
267	+ clj(1,j)*dfield(ij,1) &
268	+ clj(2,j)*dfield(ij,2) &
269	+ clj(3,j)*dfield(ij,3)
270
271	END DO
272	CASE(4)
273	!$OMP DO SCHEDULE(static)
274	DO ij=1,sz
275	field(ij) = field(ij) &
276	+ clj(1,j)*dfield(ij,1) &
277	+ clj(2,j)*dfield(ij,2) &
278	+ clj(3,j)*dfield(ij,3) &
279	+ clj(4,j)*dfield(ij,4)
280	END DO
281	CASE(5)
282	!$OMP DO SCHEDULE(static)
283	DO ij=1,sz
284	field(ij) = field(ij) &
285	+ clj(1,j)*dfield(ij,1) &
286	+ clj(2,j)*dfield(ij,2) &
287	+ clj(3,j)*dfield(ij,3) &
288	+ clj(4,j)*dfield(ij,4) &
289	+ clj(5,j)*dfield(ij,5)
290	END DO
291	END SELECT
292	CALL exit_trace()
293	END SUBROUTINE update
294
295	!---------------------------------------------- XIOS ----------------------------------------
296
297	#ifdef CPP_USING_XIOS
298
299	SUBROUTINE setup_xios() BINDC(setup_xios)
300	! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine
301	! This is the case when calling from Python after importing mpi4py
302	INTEGER :: ierr, mpi_threading_mode
303
304	PRINT *, 'Initialize XIOS and obtain our communicator'
305	CALL xios_initialize("icosagcm",return_comm=comm_icosa)
306
307	PRINT *, 'Initialize our XIOS context'
308
309	CALL xios_context_initialize("icosagcm",comm_icosa)
310	CALL xios_get_handle("icosagcm",ctx_hdl)
311	CALL xios_set_current_context(ctx_hdl)
312	END SUBROUTINE setup_xios
313
314	SUBROUTINE call_xios_set_timestep(dt) BINDC(xios_set_timestep)
315	DBL, VALUE :: dt
316	TYPE(xios_duration) :: dtime
317	dtime%second=dt
318	CALL xios_set_timestep(dtime)
319	END SUBROUTINE call_xios_set_timestep
320
321	SUBROUTINE call_xios_update_calendar(step) BINDC(xios_update_calendar)
322	INDEX, value :: step
323	CALL xios_update_calendar(step)
324	END SUBROUTINE call_xios_update_calendar
325
326	#endif
327
328	END MODULE timestep_unstructured_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats: