source: codes/icosagcm/trunk/src/physics_interface.f90 @ 472

Last change on this file since 472 was 472, checked in by ymipsl, 8 years ago

Repairing openMP :
physic_column and physic_dcmip2016 seems OK. Only for thread assigned to domain.

YM
File size: 10.7 KB
Line 
1MODULE physics_interface_mod
2
3  USE prec
4
5  PRIVATE
6
7  TYPE t_physics_inout
8     ! Input, time-independent
9     INTEGER :: ngrid
10     REAL(rstd) :: dt_phys
11     REAL(rstd), DIMENSION(:), POINTER :: Ai, lon, lat, phis
12     ! Input, time-dependent
13     REAL(rstd), DIMENSION(:,:), POINTER :: p, pk, Temp, ulon, ulat
14     REAL(rstd), DIMENSION(:,:,:), POINTER :: q
15     ! Output arrays
16     REAL(rstd), DIMENSION(:,:), POINTER :: dTemp, dulon, dulat
17     REAL(rstd), DIMENSION(:,:,:), POINTER :: dq
18  END TYPE t_physics_inout
19
20! physics_inout is used to exchange information with physics
21! Field ngrid is initialized by physics.f90/init_physics. Its other fields
22! must be defined by XX/init_physics (where XX =  e.g. physics_dcmip.f90)
23! by either pointing to internal data of the physics package
24! or by a specific allocation
25! size : (ngrid), (ngrid,llm) except p(ngrid,llm+1), (ngrid,llm,nqtot)
26
27  TYPE(t_physics_inout), SAVE :: physics_inout
28!$OMP THREADPRIVATE(physics_inout)
29 
30! pack_info contains indices used by pack/unpack routines
31! to pack together the data of all the domains managed by the MPI process
32! It is initialized by physics.f90/init_physics
33
34  TYPE t_pack_info
35     INTEGER :: ngrid, & ! number of non-halo points in that domain
36          nseg ! number of segments (contigous parts) in that domain
37     ! size and start of each segment : ij domain index, k packed index
38     INTEGER, ALLOCATABLE :: n(:), ij(:), k(:)
39  END TYPE t_pack_info
40
41  TYPE(t_pack_info), ALLOCATABLE, SAVE :: pack_info(:)
42!$OMP THREADPRIVATE(pack_info)
43
44
45  INTERFACE pack_field
46     MODULE PROCEDURE pack_2D
47     MODULE PROCEDURE pack_3D
48     MODULE PROCEDURE pack_4D
49  END INTERFACE pack_field
50
51  INTERFACE unpack_field
52     MODULE PROCEDURE unpack_2D
53     MODULE PROCEDURE unpack_3D
54     MODULE PROCEDURE unpack_4D
55  END INTERFACE unpack_field
56
57  INTERFACE pack_domain
58     MODULE PROCEDURE pack_domain_2D
59     MODULE PROCEDURE pack_domain_3D
60     MODULE PROCEDURE pack_domain_4D
61  END INTERFACE pack_domain
62
63  INTERFACE unpack_domain
64     MODULE PROCEDURE unpack_domain_2D
65     MODULE PROCEDURE unpack_domain_3D
66     MODULE PROCEDURE unpack_domain_4D
67  END INTERFACE unpack_domain
68
69  PUBLIC :: nb_extra_physics_2D, nb_extra_physics_3D, &
70       t_physics_inout, physics_inout, &
71       t_pack_info, pack_info, init_pack_before, init_pack_after, &
72       pack_domain, pack_field, unpack_domain, unpack_field, &
73       garbage_3D
74
75CONTAINS
76 
77    SUBROUTINE init_pack_before
78    USE icosa
79    IMPLICIT NONE
80    INTEGER :: ind, offset
81
82    offset=0
83    ALLOCATE(pack_info(ndomain))
84    DO ind=1,ndomain
85       IF (.NOT. assigned_domain(ind)) CYCLE
86       CALL swap_dimensions(ind)
87       CALL swap_geometry(ind)
88       CALL count_segments(domain(ind)%own, pack_info(ind))
89       pack_info(ind)%k = pack_info(ind)%k + offset
90       offset = offset + pack_info(ind)%ngrid
91    END DO
92    physics_inout%ngrid = offset
93
94  END SUBROUTINE init_pack_before
95
96  SUBROUTINE count_segments(own, info)
97    USE icosa
98    IMPLICIT NONE
99    LOGICAL, DIMENSION(:,:) :: own
100    TYPE(t_pack_info) :: info
101
102    INTEGER, DIMENSION(jjm) :: n
103    INTEGER :: ngrid, nseg, i, j, jj, k
104    INTEGER, PARAMETER :: method=4
105    SELECT CASE(method)
106    CASE(1) ! Copy all points, including halo (works)
107       info%nseg=1
108       info%ngrid=iim*jjm
109       ALLOCATE(info%n(1))
110       ALLOCATE(info%ij(1))
111       ALLOCATE(info%k(1))
112       info%n(1)=iim*jjm
113       info%ij(1)=1
114       info%k(1)=1
115    CASE(2) ! Copy all points, including halo, one at a time (works, slow ?)
116       info%nseg=iim*jjm
117       info%ngrid=iim*jjm
118       ALLOCATE(info%n(iim*jjm))
119       ALLOCATE(info%ij(iim*jjm))
120       ALLOCATE(info%k(iim*jjm))
121       DO jj=1,iim*jjm
122          info%n(jj) =1
123          info%ij(jj)=jj
124          info%k(jj) =jj
125       END DO
126    CASE(3) ! Copy non-halo points only, one at a time (works, slow ?)
127       n=COUNT(own,1)
128       ngrid=SUM(n)
129       info%ngrid=ngrid
130       info%nseg=ngrid
131       ALLOCATE(info%n(ngrid))
132       ALLOCATE(info%ij(ngrid))
133       ALLOCATE(info%k(ngrid))
134       jj=1
135       DO j=1,jjm
136          DO i=1,iim
137             IF(own(i,j)) THEN
138                info%n(jj)=1
139                info%k(jj)=jj
140                info%ij(jj) = iim*(j-1)+i
141                jj=jj+1
142             END IF
143          END DO
144       END DO
145       
146    CASE DEFAULT ! Copy non-halo points only, as contiguous segments (works)
147       n=0
148       n=COUNT(own,1)
149       ngrid=SUM(n)
150       info%ngrid=ngrid
151       nseg=COUNT(n>0)
152       info%nseg=nseg
153       ALLOCATE(info%n(nseg))
154       ALLOCATE(info%ij(nseg))
155       ALLOCATE(info%k(nseg))
156
157       jj=1
158       k=1
159       DO j=1,jjm
160          IF(n(j)>0) THEN
161             ! find first .TRUE. value in own(:,j)
162             DO i=1,iim
163                IF(own(i,j)) THEN
164                   info%n(jj)=n(j)
165                   info%k(jj)=k
166                   info%ij(jj) = iim*(j-1)+i
167                   IF(COUNT(own(i:i+n(j)-1,j)) /= n(j)) STOP
168                   EXIT
169                END IF
170             END DO
171             k = k + n(j)
172             jj=jj+1
173          END IF
174       END DO
175
176       IF(k-1/=ngrid) THEN
177          PRINT *, 'Total number of grid points inconsistent', k-1, ngrid
178          STOP
179       END IF
180       IF(jj-1/=nseg) THEN
181          PRINT *, 'Number of segments inconsistent', jj-1, nseg
182          STOP
183       END IF
184
185    END SELECT
186   
187    PRINT *, 'count_segments', info%nseg, info%ngrid, SUM(info%n), COUNT(own), iim*jjm
188  END SUBROUTINE count_segments
189
190  SUBROUTINE init_pack_after
191    USE icosa
192    IMPLICIT NONE
193    INTEGER :: ind, offset
194    DO ind=1,ndomain
195       IF (.NOT. assigned_domain(ind)) CYCLE
196       CALL swap_dimensions(ind)
197       CALL swap_geometry(ind)
198       CALL pack_domain_2D(pack_info(ind), Ai, physics_inout%Ai)
199       CALL pack_domain_2D(pack_info(ind), lon_i, physics_inout%lon)
200       CALL pack_domain_2D(pack_info(ind), lat_i, physics_inout%lat)
201    END DO
202  END SUBROUTINE init_pack_after
203
204!-------------------------------- Pack / Unpack 2D ---------------------------
205
206  SUBROUTINE pack_2D(f_2D, packed)
207    USE icosa
208    IMPLICIT NONE
209    TYPE(t_field),POINTER :: f_2D(:)
210    REAL(rstd)            :: packed(:)
211    REAL(rstd), POINTER   :: loc(:)
212    INTEGER :: ind
213    DO ind=1,ndomain
214       IF (.NOT. assigned_domain(ind)) CYCLE
215       loc = f_2D(ind)
216       CALL pack_domain_2D(pack_info(ind), loc, packed)
217    END DO
218  END SUBROUTINE pack_2D
219
220  SUBROUTINE unpack_2D(f_2D, packed) 
221    USE icosa
222    IMPLICIT NONE
223    TYPE(t_field),POINTER :: f_2D(:)
224    REAL(rstd)            :: packed(:)
225    REAL(rstd), POINTER   :: loc(:)
226    INTEGER :: ind
227    DO ind=1,ndomain
228       IF (.NOT. assigned_domain(ind)) CYCLE
229       loc = f_2D(ind)
230       CALL unpack_domain_2D(pack_info(ind), loc, packed)
231    END DO
232  END SUBROUTINE unpack_2D
233
234  SUBROUTINE pack_domain_2D(info, loc, glob)
235    USE icosa
236    IMPLICIT NONE
237    TYPE(t_pack_info) :: info
238    REAL(rstd), DIMENSION(:) :: loc, glob
239    INTEGER :: jj,n,k,ij
240    DO jj=1, info%nseg
241       n = info%n(jj)-1
242       ij = info%ij(jj)
243       k = info%k(jj)
244       glob(k:k+n) = loc(ij:ij+n)
245    END DO
246  END SUBROUTINE pack_domain_2D
247
248  SUBROUTINE unpack_domain_2D(info, loc, glob)
249    IMPLICIT NONE
250    TYPE(t_pack_info) :: info
251    REAL(rstd), DIMENSION(:) :: loc, glob
252    INTEGER :: jj,n,k,ij
253    DO jj=1, info%nseg
254       n = info%n(jj)-1
255       ij = info%ij(jj)
256       k = info%k(jj)
257       loc(ij:ij+n) = glob(k:k+n)
258    END DO
259  END SUBROUTINE unpack_domain_2D
260
261!-------------------------------- Pack / Unpack 3D ---------------------------
262
263  SUBROUTINE pack_3D(f_3D, packed)
264    USE icosa
265    IMPLICIT NONE
266    TYPE(t_field),POINTER :: f_3D(:)
267    REAL(rstd)            :: packed(:,:)
268    REAL(rstd), POINTER   :: loc(:,:)
269    INTEGER :: ind
270    DO ind=1,ndomain
271       IF (.NOT. assigned_domain(ind)) CYCLE
272       loc = f_3D(ind)
273       CALL pack_domain_3D(pack_info(ind), loc, packed)
274    END DO
275  END SUBROUTINE pack_3D
276
277  SUBROUTINE unpack_3D(f_3D, packed) 
278    USE icosa
279    IMPLICIT NONE
280    TYPE(t_field),POINTER :: f_3D(:)
281    REAL(rstd)            :: packed(:,:)
282    REAL(rstd), POINTER   :: loc(:,:)
283    INTEGER :: ind
284    DO ind=1,ndomain
285       IF (.NOT. assigned_domain(ind)) CYCLE
286       loc = f_3D(ind)
287       CALL unpack_domain_3D(pack_info(ind), loc, packed)
288    END DO
289  END SUBROUTINE unpack_3D
290
291  SUBROUTINE pack_domain_3D(info, loc, glob)
292    IMPLICIT NONE
293    TYPE(t_pack_info) :: info
294    REAL(rstd), DIMENSION(:,:) :: loc, glob
295    INTEGER :: jj,n,k,ij
296    DO jj=1, info%nseg
297       n = info%n(jj)-1
298       ij = info%ij(jj)
299       k = info%k(jj)
300       glob(k:k+n,:) = loc(ij:ij+n,:)
301    END DO
302  END SUBROUTINE pack_domain_3D
303
304  SUBROUTINE unpack_domain_3D(info, loc, glob)
305    IMPLICIT NONE
306    TYPE(t_pack_info) :: info
307    REAL(rstd), DIMENSION(:,:) :: loc, glob
308    INTEGER :: jj,n,k,ij
309    DO jj=1, info%nseg
310       n = info%n(jj)-1
311       ij = info%ij(jj)
312       k = info%k(jj)
313       loc(ij:ij+n,:) = glob(k:k+n,:)
314    END DO   
315  END SUBROUTINE unpack_domain_3D
316
317  SUBROUTINE garbage_3D(loc,own)
318    USE icosa
319    IMPLICIT NONE
320    LOGICAL :: own(iim,jjm)
321    REAL(rstd) :: loc(iim*jjm,llm)
322    INTEGER :: i,j,ij
323    ! write garbage in non-owned points
324    DO j=1,jjm
325       DO i=1,iim
326          IF(.NOT.own(i,j)) THEN
327             ij=iim*(j-1)+i
328             loc(ij,:)=-1e30
329          END IF
330       END DO
331    END DO   
332  END SUBROUTINE garbage_3D
333
334!-------------------------------- Pack / Unpack 4D ---------------------------
335
336  SUBROUTINE pack_4D(f_4D, packed)
337    USE icosa
338    IMPLICIT NONE
339    TYPE(t_field),POINTER :: f_4D(:)
340    REAL(rstd)            :: packed(:,:,:)
341    REAL(rstd), POINTER   :: loc(:,:,:)
342    INTEGER :: ind
343    DO ind=1,ndomain
344       IF (.NOT. assigned_domain(ind)) CYCLE
345       loc = f_4D(ind)
346       CALL pack_domain_4D(pack_info(ind), loc, packed)
347    END DO
348  END SUBROUTINE pack_4D
349
350  SUBROUTINE unpack_4D(f_4D, packed) 
351    USE icosa
352    IMPLICIT NONE
353    TYPE(t_field),POINTER :: f_4D(:)
354    REAL(rstd)            :: packed(:,:,:)
355    REAL(rstd), POINTER   :: loc(:,:,:)
356    INTEGER :: ind
357    DO ind=1,ndomain
358       IF (.NOT. assigned_domain(ind)) CYCLE
359       loc = f_4D(ind)
360       CALL unpack_domain_4D(pack_info(ind), loc, packed)
361    END DO
362  END SUBROUTINE unpack_4D
363
364  SUBROUTINE pack_domain_4D(info, loc, glob)
365    IMPLICIT NONE
366    TYPE(t_pack_info) :: info
367    REAL(rstd), DIMENSION(:,:,:) :: loc, glob
368    INTEGER :: jj,n,k,ij
369    DO jj=1, info%nseg
370       n = info%n(jj)-1
371       ij = info%ij(jj)
372       k = info%k(jj)
373       glob(k:k+n,:,:) = loc(ij:ij+n,:,:)
374    END DO
375  END SUBROUTINE pack_domain_4D
376
377  SUBROUTINE unpack_domain_4D(info, loc, glob)
378    IMPLICIT NONE
379    TYPE(t_pack_info) :: info
380    REAL(rstd), DIMENSION(:,:,:) :: loc, glob
381    INTEGER :: jj,n,k,ij
382    DO jj=1, info%nseg
383       n = info%n(jj)-1
384       ij = info%ij(jj)
385       k = info%k(jj)
386       loc(ij:ij+n,:,:) = glob(k:k+n,:,:)
387    END DO
388  END SUBROUTINE unpack_domain_4D
389
390END MODULE physics_interface_mod
Note: See TracBrowser for help on using the repository browser.