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, 8*sz*(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_out*xios_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 |
---|