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

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

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