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

Last change on this file since 214 was 214, checked in by dubos, 10 years ago

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