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

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

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