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

Last change on this file since 784 was 784, checked in by jisesh, 6 years ago

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