source: codes/icosagcm/devel/src/physics/physics_interface.f90 @ 739

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

devel : small cleanup in idealized physics
File size: 11.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, ngrid
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    ngrid=offset
95    ! Input                                                                                                     
96    ALLOCATE(physics_inout%Ai(ngrid))
97    ALLOCATE(physics_inout%lon(ngrid))
98    ALLOCATE(physics_inout%lat(ngrid))
99    ALLOCATE(physics_inout%phis(ngrid))
100    ALLOCATE(physics_inout%p(ngrid,llm+1))
101    ALLOCATE(physics_inout%pk(ngrid,llm))
102    ALLOCATE(physics_inout%Temp(ngrid,llm))
103    ALLOCATE(physics_inout%ulon(ngrid,llm))
104    ALLOCATE(physics_inout%ulat(ngrid,llm))
105    ALLOCATE(physics_inout%q(ngrid,llm,nqtot))
106    ! Output (tendencies)
107    ALLOCATE(physics_inout%dTemp(ngrid,llm))
108    ALLOCATE(physics_inout%dulon(ngrid,llm))
109    ALLOCATE(physics_inout%dulat(ngrid,llm))
110    ALLOCATE(physics_inout%dq(ngrid,llm,nqtot))
111
112  END SUBROUTINE init_pack_before
113
114  SUBROUTINE count_segments(own, info)
115    USE icosa
116    IMPLICIT NONE
117    LOGICAL, DIMENSION(:,:) :: own
118    TYPE(t_pack_info) :: info
119
120    INTEGER, DIMENSION(jjm) :: n
121    INTEGER :: ngrid, nseg, i, j, jj, k
122    INTEGER, PARAMETER :: method=4
123    SELECT CASE(method)
124    CASE(1) ! Copy all points, including halo (works)
125       info%nseg=1
126       info%ngrid=iim*jjm
127       ALLOCATE(info%n(1))
128       ALLOCATE(info%ij(1))
129       ALLOCATE(info%k(1))
130       info%n(1)=iim*jjm
131       info%ij(1)=1
132       info%k(1)=1
133    CASE(2) ! Copy all points, including halo, one at a time (works, slow ?)
134       info%nseg=iim*jjm
135       info%ngrid=iim*jjm
136       ALLOCATE(info%n(iim*jjm))
137       ALLOCATE(info%ij(iim*jjm))
138       ALLOCATE(info%k(iim*jjm))
139       DO jj=1,iim*jjm
140          info%n(jj) =1
141          info%ij(jj)=jj
142          info%k(jj) =jj
143       END DO
144    CASE(3) ! Copy non-halo points only, one at a time (works, slow ?)
145       n=0
146       n(jj_begin:jj_end)=COUNT(own(ii_begin:ii_end,jj_begin:jj_end),1)
147       ngrid=SUM(n)
148       info%ngrid=ngrid
149       info%nseg=ngrid
150       ALLOCATE(info%n(ngrid))
151       ALLOCATE(info%ij(ngrid))
152       ALLOCATE(info%k(ngrid))
153       jj=1
154       DO j=1,jjm
155          DO i=1,iim
156             IF(own(i,j)) THEN
157                info%n(jj)=1
158                info%k(jj)=jj
159                info%ij(jj) = iim*(j-1)+i
160                jj=jj+1
161             END IF
162          END DO
163       END DO
164       
165    CASE DEFAULT ! Copy non-halo points only, as contiguous segments (works)
166       n=0
167       n(jj_begin:jj_end)=COUNT(own(ii_begin:ii_end,jj_begin:jj_end),1)
168       ngrid=SUM(n)
169       info%ngrid=ngrid
170       nseg=COUNT(n>0)
171       info%nseg=nseg
172       ALLOCATE(info%n(nseg))
173       ALLOCATE(info%ij(nseg))
174       ALLOCATE(info%k(nseg))
175       info%n(:)=0
176       info%ij(:)=0
177       info%k(:)=0
178       
179       jj=1
180       k=1
181       DO j=jj_begin,jj_end
182          IF(n(j)>0) THEN
183             ! find first .TRUE. value in own(:,j)
184             DO i=ii_begin,ii_end
185                IF(own(i,j)) THEN
186                   info%n(jj)=n(j)
187                   info%k(jj)=k
188                   info%ij(jj) = iim*(j-1)+i
189                   IF(COUNT(own(i:i+n(j)-1,j)) /= n(j)) STOP
190                   EXIT
191                END IF
192             END DO
193             k = k + n(j)
194             jj=jj+1
195          END IF
196       END DO
197
198       IF(k-1/=ngrid) THEN
199          PRINT *, 'Total number of grid points inconsistent', k-1, ngrid
200          STOP
201       END IF
202       IF(jj-1/=nseg) THEN
203          PRINT *, 'Number of segments inconsistent', jj-1, nseg
204          STOP
205       END IF
206
207    END SELECT
208   
209    PRINT *, 'count_segments', info%nseg, info%ngrid, SUM(info%n), COUNT(own), iim*jjm
210  END SUBROUTINE count_segments
211
212  SUBROUTINE init_pack_after
213    USE icosa
214    IMPLICIT NONE
215    INTEGER :: ind, offset
216    DO ind=1,ndomain
217       IF (.NOT. assigned_domain(ind)) CYCLE
218       CALL swap_dimensions(ind)
219       CALL swap_geometry(ind)
220       CALL pack_domain_2D(pack_info(ind), Ai, physics_inout%Ai)
221       CALL pack_domain_2D(pack_info(ind), lon_i, physics_inout%lon)
222       CALL pack_domain_2D(pack_info(ind), lat_i, physics_inout%lat)
223    END DO
224  END SUBROUTINE init_pack_after
225
226!-------------------------------- Pack / Unpack 2D ---------------------------
227
228  SUBROUTINE pack_2D(f_2D, packed)
229    USE icosa
230    IMPLICIT NONE
231    TYPE(t_field),POINTER :: f_2D(:)
232    REAL(rstd)            :: packed(:)
233    REAL(rstd), POINTER   :: loc(:)
234    INTEGER :: ind
235    DO ind=1,ndomain
236       IF (.NOT. assigned_domain(ind)) CYCLE
237       loc = f_2D(ind)
238       CALL pack_domain_2D(pack_info(ind), loc, packed)
239    END DO
240  END SUBROUTINE pack_2D
241
242  SUBROUTINE unpack_2D(f_2D, packed) 
243    USE icosa
244    IMPLICIT NONE
245    TYPE(t_field),POINTER :: f_2D(:)
246    REAL(rstd)            :: packed(:)
247    REAL(rstd), POINTER   :: loc(:)
248    INTEGER :: ind
249    DO ind=1,ndomain
250       IF (.NOT. assigned_domain(ind)) CYCLE
251       loc = f_2D(ind)
252       CALL unpack_domain_2D(pack_info(ind), loc, packed)
253    END DO
254  END SUBROUTINE unpack_2D
255
256  SUBROUTINE pack_domain_2D(info, loc, glob)
257    USE icosa
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       glob(k:k+n) = loc(ij:ij+n)
267    END DO
268  END SUBROUTINE pack_domain_2D
269
270  SUBROUTINE unpack_domain_2D(info, loc, glob)
271    IMPLICIT NONE
272    TYPE(t_pack_info) :: info
273    REAL(rstd), DIMENSION(:) :: loc, glob
274    INTEGER :: jj,n,k,ij
275    DO jj=1, info%nseg
276       n = info%n(jj)-1
277       ij = info%ij(jj)
278       k = info%k(jj)
279       loc(ij:ij+n) = glob(k:k+n)
280    END DO
281  END SUBROUTINE unpack_domain_2D
282
283!-------------------------------- Pack / Unpack 3D ---------------------------
284
285  SUBROUTINE pack_3D(f_3D, packed)
286    USE icosa
287    IMPLICIT NONE
288    TYPE(t_field),POINTER :: f_3D(:)
289    REAL(rstd)            :: packed(:,:)
290    REAL(rstd), POINTER   :: loc(:,:)
291    INTEGER :: ind
292    DO ind=1,ndomain
293       IF (.NOT. assigned_domain(ind)) CYCLE
294       loc = f_3D(ind)
295       CALL pack_domain_3D(pack_info(ind), loc, packed)
296    END DO
297  END SUBROUTINE pack_3D
298
299  SUBROUTINE unpack_3D(f_3D, packed) 
300    USE icosa
301    IMPLICIT NONE
302    TYPE(t_field),POINTER :: f_3D(:)
303    REAL(rstd)            :: packed(:,:)
304    REAL(rstd), POINTER   :: loc(:,:)
305    INTEGER :: ind
306    DO ind=1,ndomain
307       IF (.NOT. assigned_domain(ind)) CYCLE
308       loc = f_3D(ind)
309       CALL unpack_domain_3D(pack_info(ind), loc, packed)
310    END DO
311  END SUBROUTINE unpack_3D
312
313  SUBROUTINE pack_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       glob(k:k+n,:) = loc(ij:ij+n,:)
323    END DO
324  END SUBROUTINE pack_domain_3D
325
326  SUBROUTINE unpack_domain_3D(info, loc, glob)
327    IMPLICIT NONE
328    TYPE(t_pack_info) :: info
329    REAL(rstd), DIMENSION(:,:) :: loc, glob
330    INTEGER :: jj,n,k,ij
331    DO jj=1, info%nseg
332       n = info%n(jj)-1
333       ij = info%ij(jj)
334       k = info%k(jj)
335       loc(ij:ij+n,:) = glob(k:k+n,:)
336    END DO   
337  END SUBROUTINE unpack_domain_3D
338
339  SUBROUTINE garbage_3D(loc,own)
340    USE icosa
341    IMPLICIT NONE
342    LOGICAL :: own(iim,jjm)
343    REAL(rstd) :: loc(iim*jjm,llm)
344    INTEGER :: i,j,ij
345    ! write garbage in non-owned points
346    DO j=1,jjm
347       DO i=1,iim
348          IF(.NOT.own(i,j)) THEN
349             ij=iim*(j-1)+i
350             loc(ij,:)=-1e30
351          END IF
352       END DO
353    END DO   
354  END SUBROUTINE garbage_3D
355
356!-------------------------------- Pack / Unpack 4D ---------------------------
357
358  SUBROUTINE pack_4D(f_4D, packed)
359    USE icosa
360    IMPLICIT NONE
361    TYPE(t_field),POINTER :: f_4D(:)
362    REAL(rstd)            :: packed(:,:,:)
363    REAL(rstd), POINTER   :: loc(:,:,:)
364    INTEGER :: ind
365    DO ind=1,ndomain
366       IF (.NOT. assigned_domain(ind)) CYCLE
367       loc = f_4D(ind)
368       CALL pack_domain_4D(pack_info(ind), loc, packed)
369    END DO
370  END SUBROUTINE pack_4D
371
372  SUBROUTINE unpack_4D(f_4D, packed) 
373    USE icosa
374    IMPLICIT NONE
375    TYPE(t_field),POINTER :: f_4D(:)
376    REAL(rstd)            :: packed(:,:,:)
377    REAL(rstd), POINTER   :: loc(:,:,:)
378    INTEGER :: ind
379    DO ind=1,ndomain
380       IF (.NOT. assigned_domain(ind)) CYCLE
381       loc = f_4D(ind)
382       CALL unpack_domain_4D(pack_info(ind), loc, packed)
383    END DO
384  END SUBROUTINE unpack_4D
385
386  SUBROUTINE pack_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       glob(k:k+n,:,:) = loc(ij:ij+n,:,:)
396    END DO
397  END SUBROUTINE pack_domain_4D
398
399  SUBROUTINE unpack_domain_4D(info, loc, glob)
400    IMPLICIT NONE
401    TYPE(t_pack_info) :: info
402    REAL(rstd), DIMENSION(:,:,:) :: loc, glob
403    INTEGER :: jj,n,k,ij
404    DO jj=1, info%nseg
405       n = info%n(jj)-1
406       ij = info%ij(jj)
407       k = info%k(jj)
408       loc(ij:ij+n,:,:) = glob(k:k+n,:,:)
409    END DO
410  END SUBROUTINE unpack_domain_4D
411
412END MODULE physics_interface_mod
Note: See TracBrowser for help on using the repository browser.