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

Last change on this file since 832 was 832, checked in by dubos, 5 years ago

devel : fix pure-MPI build
File size: 12.0 KB
Line 
1MODULE timestep_unstructured_mod
2  USE ISO_C_BINDING
3#ifdef CPP_USING_XIOS
4  USE xios
5#endif
6  USE caldyn_unstructured_mod
7  USE transfer_unstructured_mod
8  IMPLICIT NONE
9  PRIVATE
10  SAVE
11
12#ifdef CPP_USING_XIOS
13  TYPE(xios_context) :: ctx_hdl
14#endif
15
16CONTAINS
17
18#include "unstructured.h90"
19
20#define HASNAN(field) (ANY(.NOT.ABS(field)<1e20))
21
22   SUBROUTINE scalar_laplacian(s,laps) BINDC(scalar_laplacian)
23     FIELD_MASS   :: s,laps ! IN, OUT
24     FIELD_U      :: grads  ! BUF
25     !$OMP PARALLEL NUM_THREADS(nb_threads)
26     CALL update_halo(transfer_primal, s)
27     CALL compute_scalar_laplacian(s,grads,laps)
28     !$OMP END PARALLEL
29   END SUBROUTINE scalar_laplacian
30
31   SUBROUTINE curl_laplacian(u,lapu) BINDC(curl_laplacian)
32     FIELD_U   :: u,lapu ! IN, OUT
33     FIELD_Z   :: qv  ! BUF
34     !$OMP PARALLEL NUM_THREADS(nb_threads)
35     CALL update_halo(transfer_edge, u)
36     CALL compute_curl_laplacian(u,qv,lapu)
37     !$OMP END PARALLEL
38   END SUBROUTINE curl_laplacian
39
40  !
41  SUBROUTINE caldyn_unstructured(tau, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! IN : flow state
42                                 theta,ps,pk,hflux,qv, &    ! OUT : diags (except surface geopot : IN)
43                                 dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, &
44                                 dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(caldyn_unstructured) ! OUT : tendencies
45    NUM, VALUE :: tau
46    FIELD_MASS   :: rhodz, drhodz, pk, berni          ! IN, OUT, DIAG
47    FIELD_THETA  :: theta_rhodz, dtheta_rhodz, theta  ! IN, OUT, DIAG
48    FIELD_GEOPOT :: wflux, w, geopot, &               ! DIAG, INOUT
49         dPhi_fast, dPhi_slow, dW_fast, dW_slow       ! OUT
50    FIELD_U      :: u,du_fast,du_slow,hflux,qu        ! INOUT,OUT,OUT,DIAG
51    FIELD_Z      :: qv                                ! DIAG
52    FIELD_PS     :: ps,dmass_col,mass_col             ! OUT,OUT,IN (if eta_mass) or OUT,UNUSED,UNUSED (if eta_lag)
53    ! buffers for fields that need to be shared among OpenMP threads
54    ! vertical size is allocated as llm+1 even if only llm is needed
55    FIELD_GEOPOT :: bufm1, bufm2, bufm3
56    FIELD_UL     :: bufu1, bufu2, bufu3,bufu4
57
58    TIME          :: time1,time2
59    INTEGER :: ij
60   
61    time1=OMP_GET_WTIME()
62   
63    IF(hydrostatic) THEN
64       
65       !$OMP PARALLEL NUM_THREADS(nb_threads)
66       !$OMP DO SCHEDULE(STATIC)
67       DO ij=1,edge_num
68          du_fast(:,ij)=0.
69          du_slow(:,ij)=0.
70       END DO
71       !$OMP END DO
72       CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
73       CALL compute_geopot(rhodz,theta, ps,pk,geopot)
74       
75       CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
76       
77       CALL compute_pvort_only(rhodz,u,qv,qu)
78       CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow)
79       CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow)
80       IF(caldyn_eta == eta_mass) THEN
81          CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1)
82       END IF
83       !$OMP END PARALLEL
84       
85    ELSE ! NH
86       DO ij=1,edge_num
87          du_fast(:,ij)=0.
88          du_slow(:,ij)=0.
89       END DO
90       DO ij=1,primal_num
91          wflux(1,ij)=0.
92          wflux(llm+1,ij)=0.
93       END DO
94       CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
95       CALL compute_caldyn_solver(tau,rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W,dPhi_fast,dW_fast,du_fast)
96       CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
97       CALL compute_pvort_only(rhodz,u,qv,qu)
98       CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3,bufu1,bufu2,bufu3,bufu4, hflux,du_slow,dPhi_slow,dW_slow) 
99       CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow)
100       IF(caldyn_eta == eta_mass) THEN
101          CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1)
102          CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, du_slow,dPhi_slow,dW_slow)
103       END IF
104    END IF
105   
106    time2=OMP_GET_WTIME()
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    NUM2(       primal_num, max_nb_stage) :: dmass_col    ! OUT
127    NUM3(llm,   primal_num, max_nb_stage) :: drhodz   ! OUT
128    NUM4(llm,   primal_num, nqdyn, max_nb_stage) :: dtheta_rhodz ! OUT
129    NUM3(llm,   edge_num,   max_nb_stage) :: du_fast,du_slow ! OUT
130    NUM3(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    TIME       :: time1,time2
137    INTEGER :: step, stage, iq, ij
138   
139    time1=OMP_GET_WTIME()
140   
141    !$OMP PARALLEL NUM_THREADS(nb_threads)
142   
143    !$OMP DO SCHEDULE(STATIC)
144    DO ij=1,primal_num
145       wflux(1,ij)=0.
146       wflux(llm+1,ij)=0.
147    END DO
148    !$OMP END DO
149
150    DO step=1, nstep
151       DO stage=1, nb_stage
152         
153          CALL update_halo(transfer_primal, rhodz)
154          DO iq=1, nqdyn
155             CALL update_halo(transfer_primal, theta_rhodz(:,:,iq))
156          END DO
157
158          IF(hydrostatic) THEN
159             
160             !$OMP DO SCHEDULE(STATIC)
161             DO ij=1,edge_num
162                du_fast(:,ij,stage)=0.
163                du_slow(:,ij,stage)=0.
164             END DO
165             !$OMP 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    IF(time2>time1) elapsed = elapsed + time2-time1
233
234    CALL print_trace()
235   
236  END SUBROUTINE ARK_step
237 
238 
239  SUBROUTINE update(j,sz,clj,field,dfield)
240    INTEGER :: j, sz ! stage in ARK scheme, field size
241    NUM2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau
242    NUM1(sz) :: field
243    NUM2(sz, max_nb_stage) :: dfield
244    !
245    INTEGER :: l, ij
246    CALL enter_trace(id_update, 8*sz*(j+1) )
247    !  PRINT *, clj(1:j,j)
248    SELECT CASE(j)
249    CASE(1)
250       !$OMP DO SCHEDULE(static)
251       DO ij=1,sz
252          field(ij) = field(ij) &
253               + clj(1,j)*dfield(ij,1)
254       END DO
255    CASE(2)
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       END DO
262    CASE(3)
263       !$OMP DO SCHEDULE(static)
264       DO ij=1,sz
265          field(ij) = field(ij) &
266               + clj(1,j)*dfield(ij,1) &
267               + clj(2,j)*dfield(ij,2) &
268               + clj(3,j)*dfield(ij,3)
269         
270       END DO
271    CASE(4)
272       !$OMP DO SCHEDULE(static)
273       DO ij=1,sz
274          field(ij) = field(ij) &
275               + clj(1,j)*dfield(ij,1) &
276               + clj(2,j)*dfield(ij,2) &
277               + clj(3,j)*dfield(ij,3) &
278               + clj(4,j)*dfield(ij,4)
279       END DO
280    CASE(5)
281       !$OMP DO SCHEDULE(static)
282       DO ij=1,sz
283          field(ij) = field(ij) &
284               + clj(1,j)*dfield(ij,1) &
285               + clj(2,j)*dfield(ij,2) &
286               + clj(3,j)*dfield(ij,3) &
287               + clj(4,j)*dfield(ij,4) &
288               + clj(5,j)*dfield(ij,5)
289       END DO
290    END SELECT
291    CALL exit_trace()
292  END SUBROUTINE update
293 
294  !---------------------------------------------- XIOS ----------------------------------------
295
296#ifdef CPP_USING_XIOS
297 
298  SUBROUTINE setup_xios() BINDC(setup_xios)
299    ! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine
300    ! This is the case when calling from Python after importing mpi4py
301    INTEGER :: ierr, mpi_threading_mode
302    PRINT *, 'Initialize XIOS and obtain our communicator'
303    CALL xios_initialize("icosagcm",return_comm=comm_icosa)
304    CALL xios_context_initialize("icosagcm",comm_icosa)
305    CALL xios_get_handle("icosagcm",ctx_hdl)
306    CALL xios_set_current_context(ctx_hdl)
307  END SUBROUTINE setup_xios
308
309  SUBROUTINE call_xios_set_timestep(dt) BINDC(xios_set_timestep)
310    DBL, VALUE :: dt
311    TYPE(xios_duration)      :: dtime
312    dtime%second=dt
313    CALL xios_set_timestep(dtime)
314  END SUBROUTINE call_xios_set_timestep
315
316  SUBROUTINE call_xios_update_calendar(step) BINDC(xios_update_calendar)
317    INDEX, value :: step
318    CALL xios_update_calendar(step)
319  END SUBROUTINE call_xios_update_calendar
320
321#endif
322
323END MODULE timestep_unstructured_mod
Note: See TracBrowser for help on using the repository browser.