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

Last change on this file since 681 was 681, checked in by dubos, 6 years ago

devel/unstructured : local mesh setup + halo exchange
File size: 13.3 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  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
27CONTAINS
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, 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          CALL update_halo(transfer_primal, theta_rhodz(:,:,1)) ! FIXME
156
157          IF(hydrostatic) THEN
158             
159             !$OMP DO SCHEDULE(STATIC)
160             DO ij=1,edge_num
161                du_fast(:,ij,stage)=0.
162                du_slow(:,ij,stage)=0.
163             END DO
164             
165             CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
166             CALL compute_geopot(rhodz,theta, ps,pk,geopot)
167             
168             CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u)
169             CALL update_halo(transfer_edge, u)
170             CALL compute_pvort_only(rhodz,u,qv,qu)
171             CALL update_halo(transfer_edge, qu)
172
173             CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage))
174             CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz(:,:,stage), dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage))
175             IF(caldyn_eta == eta_mass) THEN
176                CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
177                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),bufu1)
178             END IF
179             
180          ELSE ! NH
181             CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
182             CALL compute_caldyn_solver(tauj(stage),rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W, &
183                  dPhi_fast(:,:,stage), dW_fast(:,:,stage), du_fast(:,:,stage))
184             CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage),u)
185             CALL compute_pvort_only(rhodz,u,qv,qu)
186             CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3,bufu1,bufu2,bufu3,bufu4, hflux, &
187                  du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage))
188             CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage))
189             IF(caldyn_eta == eta_mass) THEN
190                CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
191                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), bufu1)
192                CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, &
193                     du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) )
194             END IF
195          END IF ! NH
196         
197          ! FIXME : mass_col is computed from rhodz
198          ! so that the DOFs are the same whatever caldyn_eta
199          ! in DYNAMICO mass_col is prognosed rather than rhodz
200         
201#define UPDATE(clj,field,dfield) update(stage,SIZE(field),clj,field,dfield)
202         
203          CALL UPDATE(cslj, rhodz, drhodz)
204          CALL UPDATE(cslj, theta_rhodz, dtheta_rhodz)
205          CALL UPDATE(cslj, u, du_slow)
206          CALL UPDATE(cflj, u, du_fast)
207         
208          IF(.NOT.hydrostatic) THEN
209             CALL UPDATE(cslj, geopot, dPhi_slow)
210             CALL UPDATE(cflj, geopot, dPhi_fast)
211             CALL UPDATE(cslj, W, dW_slow)
212             CALL UPDATE(cflj, W, dW_fast)
213          END IF
214         
215       END DO
216    END DO
217    !$OMP END PARALLEL
218   
219    time2=OMP_GET_WTIME()
220    !  CALL CPU_TIME(time2)
221    IF(time2>time1) elapsed = elapsed + time2-time1
222
223    CALL print_trace()
224   
225  END SUBROUTINE ARK_step
226 
227 
228  SUBROUTINE update(j,sz,clj,field,dfield)
229    INTEGER :: j, sz ! stage in ARK scheme, field size
230    DOUBLE2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau
231    DOUBLE1(sz) :: field
232    DOUBLE2(sz, max_nb_stage) :: dfield
233    !
234    INTEGER :: l, ij
235    CALL enter_trace(id_update, 8*sz*(j+1) )
236    !  PRINT *, clj(1:j,j)
237    SELECT CASE(j)
238    CASE(1)
239       !$OMP DO SCHEDULE(static)
240       DO ij=1,sz
241          field(ij) = field(ij) &
242               + clj(1,j)*dfield(ij,1)
243       END DO
244    CASE(2)
245       !$OMP DO SCHEDULE(static)
246       DO ij=1,sz
247          field(ij) = field(ij) &
248               + clj(1,j)*dfield(ij,1) &
249               + clj(2,j)*dfield(ij,2)
250       END DO
251    CASE(3)
252       !$OMP DO SCHEDULE(static)
253       DO ij=1,sz
254          field(ij) = field(ij) &
255               + clj(1,j)*dfield(ij,1) &
256               + clj(2,j)*dfield(ij,2) &
257               + clj(3,j)*dfield(ij,3)
258         
259       END DO
260    CASE(4)
261       !$OMP DO SCHEDULE(static)
262       DO ij=1,sz
263          field(ij) = field(ij) &
264               + clj(1,j)*dfield(ij,1) &
265               + clj(2,j)*dfield(ij,2) &
266               + clj(3,j)*dfield(ij,3) &
267               + clj(4,j)*dfield(ij,4)
268       END DO
269    CASE(5)
270       !$OMP DO SCHEDULE(static)
271       DO ij=1,sz
272          field(ij) = field(ij) &
273               + clj(1,j)*dfield(ij,1) &
274               + clj(2,j)*dfield(ij,2) &
275               + clj(3,j)*dfield(ij,3) &
276               + clj(4,j)*dfield(ij,4) &
277               + clj(5,j)*dfield(ij,5)
278       END DO
279    END SELECT
280    CALL exit_trace()
281  END SUBROUTINE update
282 
283  !---------------------------------------------- XIOS ----------------------------------------
284
285#ifdef CPP_USING_XIOS
286 
287  SUBROUTINE setup_xios() BINDC(setup_xios)
288    ! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine
289    ! This is the case when calling from Python after importing mpi4py
290    INTEGER :: ierr, mpi_threading_mode
291   
292    PRINT *, 'Initialize XIOS and obtain our communicator'
293    CALL xios_initialize("icosagcm",return_comm=comm_icosa)
294   
295    PRINT *, 'Initialize our XIOS context'
296   
297    CALL xios_context_initialize("icosagcm",comm_icosa)
298    CALL xios_get_handle("icosagcm",ctx_hdl)
299    CALL xios_set_current_context(ctx_hdl)
300    !   CALL xios_set_axis_attr("lev",n_glo=llm ,value=lev_value) ;
301    !   CALL xios_set_axis_attr("levp1",n_glo=llm+1 ,value=lev_valuep1) ;
302    !   CALL xios_set_axis_attr("nq",n_glo=nqtot, value=nq_value) ;
303    !   CALL xios_set_domaingroup_attr("i",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
304    !   CALL xios_set_domaingroup_attr("i", data_dim=1, type='unstructured' , nvertex=6, i_index=ind_glo)
305    !   CALL xios_set_domaingroup_attr("i",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
306    !   CALL xios_set_domain_attr("u",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
307    !   CALL xios_set_domain_attr("u", data_dim=1, type='unstructured' , nvertex=2, i_index=ind_glo)
308    !  CALL xios_set_domain_attr("u",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
309    !   CALL xios_set_domain_attr("v",ni_glo=ncell_tot, ibegin=displ, ni=ncell)
310    !   CALL xios_set_domain_attr("v", data_dim=1, type='unstructured' , nvertex=3)
311    !   CALL xios_set_domain_attr("v",lonvalue_1d=lon, latvalue_1d=lat, bounds_lon_1d=bounds_lon, bounds_lat_1d=bounds_lat)
312    !   CALL xios_set_timestep(dtime) 
313    !   CALL xios_set_fieldgroup_attr("standard_output", freq_op=itau_out*xios_timestep, freq_offset=(itau_out-1)*xios_timestep)   
314    ! CALL xios_close_context_definition()
315   
316    ! PRINT *, 'Read data'
317    ! CALL xios_recv_field(name,field)
318   
319    ! PRINT *,'Write data'
320    ! CALL xios_send_field(name,field_tmp)
321   
322   !   PRINT *, 'Flush to disk, clean up and die'
323   !   CALL xios_context_finalize
324  END SUBROUTINE setup_xios
325!
326  SUBROUTINE call_xios_set_timestep(dt) BINDC(xios_set_timestep)
327    DBL, VALUE :: dt
328    TYPE(xios_duration)      :: dtime
329    dtime%second=dt
330    CALL xios_set_timestep(dtime)
331  END SUBROUTINE call_xios_set_timestep
332
333  SUBROUTINE call_xios_update_calendar(step) BINDC(xios_update_calendar)
334    INDEX, value :: step
335    CALL xios_update_calendar(step)
336  END SUBROUTINE call_xios_update_calendar
337
338#endif
339
340END MODULE timestep_unstructured_mod
Note: See TracBrowser for help on using the repository browser.