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

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

Add DCMIP2016 physics

(first guess)

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 
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_domain_2D(pack_info(ind), lon_i, physics_inout%lon)
197       CALL pack_domain_2D(pack_info(ind), lat_i, physics_inout%lat)
198    END DO
199  END SUBROUTINE init_pack_after
200
201!-------------------------------- Pack / Unpack 2D ---------------------------
202
203  SUBROUTINE pack_2D(f_2D, packed)
204    USE icosa
205    IMPLICIT NONE
206    TYPE(t_field),POINTER :: f_2D(:)
207    REAL(rstd)            :: packed(:)
208    REAL(rstd), POINTER   :: loc(:)
209    INTEGER :: ind
210    DO ind=1,ndomain
211       IF (.NOT. assigned_domain(ind)) CYCLE
212       loc = f_2D(ind)
213       CALL pack_domain_2D(pack_info(ind), loc, packed)
214    END DO
215  END SUBROUTINE pack_2D
216
217  SUBROUTINE unpack_2D(f_2D, packed) 
218    USE icosa
219    IMPLICIT NONE
220    TYPE(t_field),POINTER :: f_2D(:)
221    REAL(rstd)            :: packed(:)
222    REAL(rstd), POINTER   :: loc(:)
223    INTEGER :: ind
224    DO ind=1,ndomain
225       IF (.NOT. assigned_domain(ind)) CYCLE
226       loc = f_2D(ind)
227       CALL unpack_domain_2D(pack_info(ind), loc, packed)
228    END DO
229  END SUBROUTINE unpack_2D
230
231  SUBROUTINE pack_domain_2D(info, loc, glob)
232    USE icosa
233    IMPLICIT NONE
234    TYPE(t_pack_info) :: info
235    REAL(rstd), DIMENSION(:) :: loc, glob
236    INTEGER :: jj,n,k,ij
237    DO jj=1, info%nseg
238       n = info%n(jj)-1
239       ij = info%ij(jj)
240       k = info%k(jj)
241       glob(k:k+n) = loc(ij:ij+n)
242    END DO
243  END SUBROUTINE pack_domain_2D
244
245  SUBROUTINE unpack_domain_2D(info, loc, glob)
246    IMPLICIT NONE
247    TYPE(t_pack_info) :: info
248    REAL(rstd), DIMENSION(:) :: loc, glob
249    INTEGER :: jj,n,k,ij
250    DO jj=1, info%nseg
251       n = info%n(jj)-1
252       ij = info%ij(jj)
253       k = info%k(jj)
254       loc(ij:ij+n) = glob(k:k+n)
255    END DO
256  END SUBROUTINE unpack_domain_2D
257
258!-------------------------------- Pack / Unpack 3D ---------------------------
259
260  SUBROUTINE pack_3D(f_3D, packed)
261    USE icosa
262    IMPLICIT NONE
263    TYPE(t_field),POINTER :: f_3D(:)
264    REAL(rstd)            :: packed(:,:)
265    REAL(rstd), POINTER   :: loc(:,:)
266    INTEGER :: ind
267    DO ind=1,ndomain
268       IF (.NOT. assigned_domain(ind)) CYCLE
269       loc = f_3D(ind)
270       CALL pack_domain_3D(pack_info(ind), loc, packed)
271    END DO
272  END SUBROUTINE pack_3D
273
274  SUBROUTINE unpack_3D(f_3D, packed) 
275    USE icosa
276    IMPLICIT NONE
277    TYPE(t_field),POINTER :: f_3D(:)
278    REAL(rstd)            :: packed(:,:)
279    REAL(rstd), POINTER   :: loc(:,:)
280    INTEGER :: ind
281    DO ind=1,ndomain
282       IF (.NOT. assigned_domain(ind)) CYCLE
283       loc = f_3D(ind)
284       CALL unpack_domain_3D(pack_info(ind), loc, packed)
285    END DO
286  END SUBROUTINE unpack_3D
287
288  SUBROUTINE pack_domain_3D(info, loc, glob)
289    IMPLICIT NONE
290    TYPE(t_pack_info) :: info
291    REAL(rstd), DIMENSION(:,:) :: loc, glob
292    INTEGER :: jj,n,k,ij
293    DO jj=1, info%nseg
294       n = info%n(jj)-1
295       ij = info%ij(jj)
296       k = info%k(jj)
297       glob(k:k+n,:) = loc(ij:ij+n,:)
298    END DO
299  END SUBROUTINE pack_domain_3D
300
301  SUBROUTINE unpack_domain_3D(info, loc, glob)
302    IMPLICIT NONE
303    TYPE(t_pack_info) :: info
304    REAL(rstd), DIMENSION(:,:) :: loc, glob
305    INTEGER :: jj,n,k,ij
306    DO jj=1, info%nseg
307       n = info%n(jj)-1
308       ij = info%ij(jj)
309       k = info%k(jj)
310       loc(ij:ij+n,:) = glob(k:k+n,:)
311    END DO   
312  END SUBROUTINE unpack_domain_3D
313
314  SUBROUTINE garbage_3D(loc,own)
315    USE icosa
316    IMPLICIT NONE
317    LOGICAL :: own(iim,jjm)
318    REAL(rstd) :: loc(iim*jjm,llm)
319    INTEGER :: i,j,ij
320    ! write garbage in non-owned points
321    DO j=1,jjm
322       DO i=1,iim
323          IF(.NOT.own(i,j)) THEN
324             ij=iim*(j-1)+i
325             loc(ij,:)=-1e30
326          END IF
327       END DO
328    END DO   
329  END SUBROUTINE garbage_3D
330
331!-------------------------------- Pack / Unpack 4D ---------------------------
332
333  SUBROUTINE pack_4D(f_4D, packed)
334    USE icosa
335    IMPLICIT NONE
336    TYPE(t_field),POINTER :: f_4D(:)
337    REAL(rstd)            :: packed(:,:,:)
338    REAL(rstd), POINTER   :: loc(:,:,:)
339    INTEGER :: ind
340    DO ind=1,ndomain
341       IF (.NOT. assigned_domain(ind)) CYCLE
342       loc = f_4D(ind)
343       CALL pack_domain_4D(pack_info(ind), loc, packed)
344    END DO
345  END SUBROUTINE pack_4D
346
347  SUBROUTINE unpack_4D(f_4D, packed) 
348    USE icosa
349    IMPLICIT NONE
350    TYPE(t_field),POINTER :: f_4D(:)
351    REAL(rstd)            :: packed(:,:,:)
352    REAL(rstd), POINTER   :: loc(:,:,:)
353    INTEGER :: ind
354    DO ind=1,ndomain
355       IF (.NOT. assigned_domain(ind)) CYCLE
356       loc = f_4D(ind)
357       CALL unpack_domain_4D(pack_info(ind), loc, packed)
358    END DO
359  END SUBROUTINE unpack_4D
360
361  SUBROUTINE pack_domain_4D(info, loc, glob)
362    IMPLICIT NONE
363    TYPE(t_pack_info) :: info
364    REAL(rstd), DIMENSION(:,:,:) :: loc, glob
365    INTEGER :: jj,n,k,ij
366    DO jj=1, info%nseg
367       n = info%n(jj)-1
368       ij = info%ij(jj)
369       k = info%k(jj)
370       glob(k:k+n,:,:) = loc(ij:ij+n,:,:)
371    END DO
372  END SUBROUTINE pack_domain_4D
373
374  SUBROUTINE unpack_domain_4D(info, loc, glob)
375    IMPLICIT NONE
376    TYPE(t_pack_info) :: info
377    REAL(rstd), DIMENSION(:,:,:) :: loc, glob
378    INTEGER :: jj,n,k,ij
379    DO jj=1, info%nseg
380       n = info%n(jj)-1
381       ij = info%ij(jj)
382       k = info%k(jj)
383       loc(ij:ij+n,:,:) = glob(k:k+n,:,:)
384    END DO
385  END SUBROUTINE unpack_domain_4D
386
387END MODULE physics_interface_mod
Note: See TracBrowser for help on using the repository browser.