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

Last change on this file since 638 was 638, checked in by dubos, 7 years ago

devel/unstructured : code reorganization
File size: 11.9 KB
Line 
1MODULE 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
26CONTAINS
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
38SUBROUTINE 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,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  DOUBLE2(llm+1, edge_num) :: wwuu
51  DBL          :: time1,time2
52  INTEGER :: ij
53
54!  CALL CPU_TIME(time1)
55  time1=OMP_GET_WTIME()
56
57  IF(hydrostatic) THEN
58
59    !$OMP PARALLEL NUM_THREADS(nb_threads)
60    !$OMP DO SCHEDULE(STATIC)
61    DO ij=1,edge_num
62      du_fast(:,ij)=0.
63      du_slow(:,ij)=0.
64    END DO
65    !$OMP END DO
66    CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
67    CALL compute_geopot(rhodz,theta, ps,pk,geopot)
68
69    CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
70
71    CALL compute_pvort_only(rhodz,u,qv,qu)
72    CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow)
73    CALL compute_coriolis(hflux,theta,qu, drhodz,dtheta_rhodz,du_slow)
74    IF(caldyn_eta == eta_mass) THEN
75       CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,wwuu)
76    END IF
77!$OMP END PARALLEL
78
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, 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
104END SUBROUTINE caldyn_unstructured
105!
106!----------------------------- Time stepping -------------------------------
107
108!
109  SUBROUTINE ARK_step(mass_col,rhodz,theta_rhodz,u,geopot,w, & ! INOUT : flow state
110       theta,ps,pk,hflux,qv, &    ! OUT : diags (except surface geopot : IN)
111       dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, &
112       dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(ARK_step) ! OUT : tendencies
113    FIELD_PS     :: mass_col, ps                     ! OUT,IN (if eta_mass) or OUT,UNUSED (if eta_lag)
114    FIELD_MASS   :: rhodz, pk, berni                 ! IN, DIAG
115    FIELD_THETA  :: theta_rhodz, theta               ! IN, DIAG
116    FIELD_U      :: u,hflux,qu                       ! INOUT,DIAG
117    FIELD_GEOPOT :: geopot, w, wflux, w_il           ! IN, INOUT, DIAG, DIAG
118    FIELD_Z      :: qv                               ! DIAG
119    DOUBLE2(       primal_num, max_nb_stage) :: dmass_col    ! OUT
120    DOUBLE3(llm,   primal_num, max_nb_stage) :: drhodz   ! OUT
121    DOUBLE4(llm,   primal_num, nqdyn, max_nb_stage) :: dtheta_rhodz ! OUT
122    DOUBLE3(llm,   edge_num,   max_nb_stage) :: du_fast,du_slow ! OUT
123    DOUBLE3(llm+1, primal_num, max_nb_stage) :: &
124         dPhi_fast, dPhi_slow, dW_fast, dW_slow      ! OUT
125    FIELD_UL     :: DePhil, v_el, G_el               ! DIAG*3
126    DOUBLE2(llm+1, edge_num) :: wwuu
127    DBL       :: time1,time2
128    INTEGER :: stage, ij
129
130    !CALL CPU_TIME(time1)
131    time1=OMP_GET_WTIME()
132
133    !$OMP PARALLEL NUM_THREADS(nb_threads)
134   
135    DO stage=1, nb_stage
136       
137       IF(hydrostatic) THEN
138
139          !$OMP DO SCHEDULE(STATIC)
140          DO ij=1,edge_num
141             du_fast(:,ij,stage)=0.
142             du_slow(:,ij,stage)=0.
143          END DO
144
145          CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
146          CALL compute_geopot(rhodz,theta, ps,pk,geopot)
147         
148          CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u)
149         
150          CALL compute_pvort_only(rhodz,u,qv,qu)
151          CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage))
152          CALL compute_coriolis(hflux,theta,qu, drhodz(:,:,stage), dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage))
153          IF(caldyn_eta == eta_mass) THEN
154             STOP ! FIXME
155             CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
156                  dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),wwuu)
157          END IF
158         
159       ELSE ! NH
160          STOP ! FIXME
161          !$OMP DO SCHEDULE(STATIC)
162          DO ij=1,primal_num
163             wflux(1,ij)=0.
164             wflux(llm+1,ij)=0.
165          END DO
166          !$OMP END DO
167          CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
168          CALL compute_caldyn_solver(tauj(stage),rhodz,theta,pk,geopot,W, &
169               dPhi_fast(:,:,stage), dW_fast(:,:,stage), du_fast(:,:,stage))
170          CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage),u)
171          CALL compute_pvort_only(rhodz,u,qv,qu)
172          CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, hflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage))
173          CALL compute_coriolis(hflux,theta,qu, drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage))
174          IF(caldyn_eta == eta_mass) THEN
175             CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
176                  dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), wwuu)
177             CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) )
178          END IF
179       END IF ! NH
180       
181      ! FIXME : mass_col is computed from rhodz
182      ! so that the DOFs are the same whatever caldyn_eta
183      ! in DYNAMICO mass_col is prognosed rather than rhodz
184
185#define UPDATE(clj,field,dfield) update(stage,SIZE(field),clj,field,dfield)
186
187       CALL UPDATE(cslj, rhodz, drhodz)
188       CALL UPDATE(cslj, theta_rhodz, dtheta_rhodz)
189       CALL UPDATE(cslj, u, du_slow)
190       CALL UPDATE(cflj, u, du_fast)
191
192       IF(.NOT.hydrostatic) THEN
193          STOP ! FIXME
194          CALL UPDATE(cslj, geopot, dPhi_slow)
195          CALL UPDATE(cflj, geopot, dPhi_fast)
196          CALL UPDATE(cslj, W, dW_slow)
197          CALL UPDATE(cflj, W, dW_fast)
198       END IF
199
200    END DO
201!$OMP END PARALLEL
202   
203  time2=OMP_GET_WTIME()
204!  CALL CPU_TIME(time2)
205  IF(time2>time1) elapsed = elapsed + time2-time1
206
207  END SUBROUTINE ARK_step
208  !
209!
210SUBROUTINE update(j,sz,clj,field,dfield)
211  INTEGER :: j, sz ! stage in ARK scheme, field size
212  DOUBLE2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau
213  DOUBLE1(sz) :: field
214  DOUBLE2(sz, max_nb_stage) :: dfield
215  !
216  INTEGER :: l, ij
217!  PRINT *, clj(1:j,j)
218  SELECT CASE(j)
219  CASE(1)
220     !$OMP DO SCHEDULE(static)
221     DO ij=1,sz
222        field(ij) = field(ij) &
223             + clj(1,j)*dfield(ij,1)
224     END DO
225  CASE(2)
226     !$OMP DO SCHEDULE(static)
227     DO ij=1,sz
228        field(ij) = field(ij) &
229             + clj(1,j)*dfield(ij,1) &
230             + clj(2,j)*dfield(ij,2)
231     END DO
232  CASE(3)
233     !$OMP DO SCHEDULE(static)
234     DO ij=1,sz
235        field(ij) = field(ij) &
236             + clj(1,j)*dfield(ij,1) &
237             + clj(2,j)*dfield(ij,2) &
238             + clj(3,j)*dfield(ij,3)
239
240     END DO
241  CASE(4)
242     !$OMP DO SCHEDULE(static)
243     DO ij=1,sz
244        field(ij) = field(ij) &
245             + clj(1,j)*dfield(ij,1) &
246             + clj(2,j)*dfield(ij,2) &
247             + clj(3,j)*dfield(ij,3) &
248             + clj(4,j)*dfield(ij,4)
249     END DO
250  CASE(5)
251     !$OMP DO SCHEDULE(static)
252     DO ij=1,sz
253        field(ij) = field(ij) &
254             + clj(1,j)*dfield(ij,1) &
255             + clj(2,j)*dfield(ij,2) &
256             + clj(3,j)*dfield(ij,3) &
257             + clj(4,j)*dfield(ij,4) &
258             + clj(5,j)*dfield(ij,5)
259     END DO
260  END SELECT
261END SUBROUTINE update
262
263!---------------------------------------------- XIOS ----------------------------------------
264
265#ifdef CPP_USING_XIOS
266
267  SUBROUTINE setup_xios() BINDC(setup_xios)
268    ! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine
269    ! This is the case when calling from Python after importing mpi4py
270    INTEGER :: ierr, mpi_threading_mode
271   
272    PRINT *, 'Initialize XIOS and obtain our communicator'
273    CALL xios_initialize("icosagcm",return_comm=comm_icosa)
274   
275    PRINT *, 'Initialize our XIOS context'
276   
277    CALL xios_context_initialize("icosagcm",comm_icosa)
278    CALL xios_get_handle("icosagcm",ctx_hdl)
279    CALL xios_set_current_context(ctx_hdl)
280    !   CALL xios_set_axis_attr("lev",n_glo=llm ,value=lev_value) ;
281    !   CALL xios_set_axis_attr("levp1",n_glo=llm+1 ,value=lev_valuep1) ;
282    !   CALL xios_set_axis_attr("nq",n_glo=nqtot, value=nq_value) ;
283    !   CALL xios_set_domaingroup_attr("i",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
284    !   CALL xios_set_domaingroup_attr("i", data_dim=1, type='unstructured' , nvertex=6, i_index=ind_glo)
285    !   CALL xios_set_domaingroup_attr("i",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
286    !   CALL xios_set_domain_attr("u",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
287    !   CALL xios_set_domain_attr("u", data_dim=1, type='unstructured' , nvertex=2, i_index=ind_glo)
288    !  CALL xios_set_domain_attr("u",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
289    !   CALL xios_set_domain_attr("v",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
290    !   CALL xios_set_domain_attr("v", data_dim=1, type='unstructured' , nvertex=3)
291    !   CALL xios_set_domain_attr("v",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
292    !   CALL xios_set_timestep(dtime) 
293    !   CALL xios_set_fieldgroup_attr("standard_output", freq_op=itau_out*xios_timestep, freq_offset=(itau_out-1)*xios_timestep)   
294    ! CALL xios_close_context_definition()
295   
296    ! PRINT *, 'Read data'
297    ! CALL xios_recv_field(name,field)
298   
299    ! PRINT *,'Write data'
300    ! CALL xios_send_field(name,field_tmp)
301   
302   !   PRINT *, 'Flush to disk, clean up and die'
303   !   CALL xios_context_finalize
304  END SUBROUTINE setup_xios
305!
306  SUBROUTINE call_xios_set_timestep(dt) BINDC(xios_set_timestep)
307    DBL, VALUE :: dt
308    TYPE(xios_duration)      :: dtime
309    dtime%second=dt
310    CALL xios_set_timestep(dtime)
311  END SUBROUTINE call_xios_set_timestep
312
313  SUBROUTINE call_xios_update_calendar(step) BINDC(xios_update_calendar)
314    INDEX, value :: step
315    CALL xios_update_calendar(step)
316  END SUBROUTINE call_xios_update_calendar
317
318#endif
319
320END MODULE timestep_unstructured_mod
Note: See TracBrowser for help on using the repository browser.