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

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

devel/unstructured : fix missing OMP END DO
File size: 11.5 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#ifdef CPP_USING_XIOS
14  TYPE(xios_context) :: ctx_hdl
15#endif
16
17CONTAINS
18
19#include "unstructured.h90"
20
21#define HASNAN(field) (ANY(.NOT.ABS(field)<1e20))
22
23  SUBROUTINE caldyn_unstructured(tau, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! IN : flow state
24                                 theta,ps,pk,hflux,qv, &    ! OUT : diags (except surface geopot : IN)
25                                 dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, &
26                                 dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(caldyn_unstructured) ! OUT : tendencies
27    NUM, VALUE :: tau
28    FIELD_MASS   :: rhodz, drhodz, pk, berni          ! IN, OUT, DIAG
29    FIELD_THETA  :: theta_rhodz, dtheta_rhodz, theta  ! IN, OUT, DIAG
30    FIELD_GEOPOT :: wflux, w, geopot, &               ! DIAG, INOUT
31         dPhi_fast, dPhi_slow, dW_fast, dW_slow       ! OUT
32    FIELD_U      :: u,du_fast,du_slow,hflux,qu        ! INOUT,OUT,OUT,DIAG
33    FIELD_Z      :: qv                                ! DIAG
34    FIELD_PS     :: ps,dmass_col,mass_col             ! OUT,OUT,IN (if eta_mass) or OUT,UNUSED,UNUSED (if eta_lag)
35    ! buffers for fields that need to be shared among OpenMP threads
36    ! vertical size is allocated as llm+1 even if only llm is needed
37    FIELD_GEOPOT :: bufm1, bufm2, bufm3
38    FIELD_UL     :: bufu1, bufu2, bufu3,bufu4
39
40    TIME          :: time1,time2
41    INTEGER :: ij
42   
43    !  CALL CPU_TIME(time1)
44    time1=OMP_GET_WTIME()
45   
46    IF(hydrostatic) THEN
47       
48       !$OMP PARALLEL NUM_THREADS(nb_threads)
49       !$OMP DO SCHEDULE(STATIC)
50       DO ij=1,edge_num
51          du_fast(:,ij)=0.
52          du_slow(:,ij)=0.
53       END DO
54       !$OMP END DO
55       CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
56       CALL compute_geopot(rhodz,theta, ps,pk,geopot)
57       
58       CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
59       
60       CALL compute_pvort_only(rhodz,u,qv,qu)
61       CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow)
62       CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow)
63       IF(caldyn_eta == eta_mass) THEN
64          CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1)
65       END IF
66       !$OMP END PARALLEL
67       
68    ELSE ! NH
69       DO ij=1,edge_num
70          du_fast(:,ij)=0.
71          du_slow(:,ij)=0.
72       END DO
73       DO ij=1,primal_num
74          wflux(1,ij)=0.
75          wflux(llm+1,ij)=0.
76       END DO
77       CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
78       CALL compute_caldyn_solver(tau,rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W,dPhi_fast,dW_fast,du_fast)
79       CALL compute_caldyn_fast(tau, pk,berni,theta,geopot, du_fast,u)
80       CALL compute_pvort_only(rhodz,u,qv,qu)
81       CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3,bufu1,bufu2,bufu3,bufu4, hflux,du_slow,dPhi_slow,dW_slow) 
82       CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz,dtheta_rhodz,du_slow)
83       IF(caldyn_eta == eta_mass) THEN
84          CALL caldyn_vert(drhodz,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du_slow,bufu1)
85          CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, du_slow,dPhi_slow,dW_slow)
86       END IF
87    END IF
88   
89    time2=OMP_GET_WTIME()
90    !  CALL CPU_TIME(time2)
91    IF(time2>time1) elapsed = elapsed + time2-time1
92
93    CALL print_trace()
94  END SUBROUTINE caldyn_unstructured
95  !
96  !----------------------------- Time stepping -------------------------------
97 
98  !
99  SUBROUTINE ARK_step(nstep, mass_col,rhodz,theta_rhodz,u,geopot,w, & ! INOUT : flow state
100       theta,ps,pk,hflux,qv, &    ! OUT : diags (except surface geopot : IN)
101       dmass_col,drhodz,dtheta_rhodz,du_fast,du_slow, &
102       dPhi_fast, dPhi_slow, dW_fast, dW_slow) BINDC(ARK_step) ! OUT : tendencies
103    INTEGER, VALUE :: nstep ! advance by nstep time steps
104    FIELD_PS     :: mass_col, ps                     ! OUT,IN (if eta_mass) or OUT,UNUSED (if eta_lag)
105    FIELD_MASS   :: rhodz, pk, berni                 ! IN, DIAG
106    FIELD_THETA  :: theta_rhodz, theta               ! IN, DIAG
107    FIELD_U      :: u,hflux,qu                       ! INOUT,DIAG*2
108    FIELD_GEOPOT :: geopot, w, wflux                 ! IN, INOUT, DIAG
109    FIELD_Z      :: qv                               ! DIAG
110    NUM2(       primal_num, max_nb_stage) :: dmass_col    ! OUT
111    NUM3(llm,   primal_num, max_nb_stage) :: drhodz   ! OUT
112    NUM4(llm,   primal_num, nqdyn, max_nb_stage) :: dtheta_rhodz ! OUT
113    NUM3(llm,   edge_num,   max_nb_stage) :: du_fast,du_slow ! OUT
114    NUM3(llm+1, primal_num, max_nb_stage) :: &
115         dPhi_fast, dPhi_slow, dW_fast, dW_slow      ! OUT
116    ! buffers for fields that need to be shared among OpenMP threads
117    ! vertical size is allocated as llm+1 even if only llm is needed
118    FIELD_GEOPOT :: bufm1, bufm2, bufm3
119    FIELD_UL     :: bufu1, bufu2, bufu3,bufu4
120    TIME       :: time1,time2
121    INTEGER :: step, stage, iq, ij
122   
123    !CALL CPU_TIME(time1)
124    time1=OMP_GET_WTIME()
125   
126    !$OMP PARALLEL NUM_THREADS(nb_threads)
127   
128    !$OMP DO SCHEDULE(STATIC)
129    DO ij=1,primal_num
130       wflux(1,ij)=0.
131       wflux(llm+1,ij)=0.
132    END DO
133    !$OMP END DO
134
135    DO step=1, nstep
136       DO stage=1, nb_stage
137         
138          CALL update_halo(transfer_primal, rhodz)
139          DO iq=1, nqdyn
140             CALL update_halo(transfer_primal, theta_rhodz(:,:,iq))
141          END DO
142
143          IF(hydrostatic) THEN
144             
145             !$OMP DO SCHEDULE(STATIC)
146             DO ij=1,edge_num
147                du_fast(:,ij,stage)=0.
148                du_slow(:,ij,stage)=0.
149             END DO
150             !$OMP END DO
151             
152             CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
153             CALL compute_geopot(rhodz,theta, ps,pk,geopot)
154             
155             CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage), u)
156             CALL update_halo(transfer_edge, u)
157             CALL compute_pvort_only(rhodz,u,qv,qu)
158             CALL update_halo(transfer_edge, qu)
159           
160             CALL compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du_slow(:,:,stage))
161             CALL compute_coriolis(hflux,theta,qu, bufu1, drhodz(:,:,stage), &
162                  dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage))
163             IF(caldyn_eta == eta_mass) THEN
164                CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
165                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage),du_slow(:,:,stage),bufu1)
166             END IF
167             
168          ELSE ! NH
169             CALL update_halo(transfer_primal, geopot)
170             CALL update_halo(transfer_primal, W)
171             
172             CALL compute_theta(mass_col,rhodz,theta_rhodz, theta)
173             CALL compute_caldyn_solver(tauj(stage),rhodz,theta, bufm1,bufm2,bufm3, pk,geopot,W, &
174                  dPhi_fast(:,:,stage), dW_fast(:,:,stage), du_fast(:,:,stage))
175
176             CALL compute_caldyn_fast(tauj(stage), pk,berni,theta,geopot, du_fast(:,:,stage),u)
177             CALL update_halo(transfer_edge, u)
178             CALL compute_pvort_only(rhodz,u,qv,qu)
179             CALL update_halo(transfer_edge, qu)
180
181             CALL compute_caldyn_slow_NH(u,rhodz,geopot,W, bufm1,bufm2,bufm3, &
182                  bufu1,bufu2,bufu3,bufu4, hflux, &
183                  du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage))
184             CALL compute_coriolis(hflux,theta,qu, bufu1, &
185                  drhodz(:,:,stage),dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage))
186             IF(caldyn_eta == eta_mass) THEN
187                CALL caldyn_vert(drhodz(:,:,stage),rhodz,theta,u, &
188                     dmass_col(:,stage),wflux,dtheta_rhodz(:,:,:,stage), du_slow(:,:,stage), bufu1)
189                CALL compute_caldyn_vert_NH(rhodz,geopot,W,wflux, bufm1,bufm2,bufm3, &
190                     du_slow(:,:,stage), dPhi_slow(:,:,stage), dW_slow(:,:,stage) )
191             END IF
192          END IF ! NH
193         
194          ! FIXME : mass_col is computed from rhodz
195          ! so that the DOFs are the same whatever caldyn_eta
196          ! in DYNAMICO mass_col is prognosed rather than rhodz
197         
198#define UPDATE(clj,field,dfield) update(stage,SIZE(field),clj,field,dfield)
199         
200          CALL UPDATE(cslj, rhodz, drhodz)
201          CALL UPDATE(cslj, theta_rhodz, dtheta_rhodz)
202          CALL UPDATE(cslj, u, du_slow)
203          CALL UPDATE(cflj, u, du_fast)
204         
205          IF(.NOT.hydrostatic) THEN
206             CALL UPDATE(cslj, geopot, dPhi_slow)
207             CALL UPDATE(cflj, geopot, dPhi_fast)
208             CALL UPDATE(cslj, W, dW_slow)
209             CALL UPDATE(cflj, W, dW_fast)
210          END IF
211         
212       END DO
213    END DO
214    !$OMP END PARALLEL
215   
216    time2=OMP_GET_WTIME()
217    !  CALL CPU_TIME(time2)
218    IF(time2>time1) elapsed = elapsed + time2-time1
219
220    CALL print_trace()
221   
222  END SUBROUTINE ARK_step
223 
224 
225  SUBROUTINE update(j,sz,clj,field,dfield)
226    INTEGER :: j, sz ! stage in ARK scheme, field size
227    NUM2(max_nb_stage,max_nb_stage) :: clj ! modified Butcher tableau
228    NUM1(sz) :: field
229    NUM2(sz, max_nb_stage) :: dfield
230    !
231    INTEGER :: l, ij
232    CALL enter_trace(id_update, 8*sz*(j+1) )
233    !  PRINT *, clj(1:j,j)
234    SELECT CASE(j)
235    CASE(1)
236       !$OMP DO SCHEDULE(static)
237       DO ij=1,sz
238          field(ij) = field(ij) &
239               + clj(1,j)*dfield(ij,1)
240       END DO
241    CASE(2)
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       END DO
248    CASE(3)
249       !$OMP DO SCHEDULE(static)
250       DO ij=1,sz
251          field(ij) = field(ij) &
252               + clj(1,j)*dfield(ij,1) &
253               + clj(2,j)*dfield(ij,2) &
254               + clj(3,j)*dfield(ij,3)
255         
256       END DO
257    CASE(4)
258       !$OMP DO SCHEDULE(static)
259       DO ij=1,sz
260          field(ij) = field(ij) &
261               + clj(1,j)*dfield(ij,1) &
262               + clj(2,j)*dfield(ij,2) &
263               + clj(3,j)*dfield(ij,3) &
264               + clj(4,j)*dfield(ij,4)
265       END DO
266    CASE(5)
267       !$OMP DO SCHEDULE(static)
268       DO ij=1,sz
269          field(ij) = field(ij) &
270               + clj(1,j)*dfield(ij,1) &
271               + clj(2,j)*dfield(ij,2) &
272               + clj(3,j)*dfield(ij,3) &
273               + clj(4,j)*dfield(ij,4) &
274               + clj(5,j)*dfield(ij,5)
275       END DO
276    END SELECT
277    CALL exit_trace()
278  END SUBROUTINE update
279 
280  !---------------------------------------------- XIOS ----------------------------------------
281
282#ifdef CPP_USING_XIOS
283 
284  SUBROUTINE setup_xios() BINDC(setup_xios)
285    ! MPI_INIT / MPI_finalize are assumed to be called BEFORE/AFTER this routine
286    ! This is the case when calling from Python after importing mpi4py
287    INTEGER :: ierr, mpi_threading_mode
288   
289    PRINT *, 'Initialize XIOS and obtain our communicator'
290    CALL xios_initialize("icosagcm",return_comm=comm_icosa)
291   
292    PRINT *, 'Initialize our XIOS context'
293   
294    CALL xios_context_initialize("icosagcm",comm_icosa)
295    CALL xios_get_handle("icosagcm",ctx_hdl)
296    CALL xios_set_current_context(ctx_hdl)
297  END SUBROUTINE setup_xios
298
299  SUBROUTINE call_xios_set_timestep(dt) BINDC(xios_set_timestep)
300    DBL, VALUE :: dt
301    TYPE(xios_duration)      :: dtime
302    dtime%second=dt
303    CALL xios_set_timestep(dtime)
304  END SUBROUTINE call_xios_set_timestep
305
306  SUBROUTINE call_xios_update_calendar(step) BINDC(xios_update_calendar)
307    INDEX, value :: step
308    CALL xios_update_calendar(step)
309  END SUBROUTINE call_xios_update_calendar
310
311#endif
312
313END MODULE timestep_unstructured_mod
Note: See TracBrowser for help on using the repository browser.