Context Navigation

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

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

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

Download in other formats: