source: codes/icosagcm/devel/src/output/set_bounds.f90 @ 1026

Last change on this file since 1026 was 993, checked in by rpennel, 5 years ago

devel : temporary hack to use both cellset and read_metric
File size: 6.8 KB
Line 
1MODULE set_bounds_mod
2  USE geometry, ONLY : swap_geometry, lon_e, lat_e
3  USE dimensions, ONLY : swap_dimensions, u_pos
4  USE math_const, ONLY : Pi
5  USE domain_mod, ONLY : t_domain, t_cellset, domloc_glo_ind
6  USE spherical_geom_mod, ONLY : xyz2lonlat
7
8  IMPLICIT NONE
9  PRIVATE
10  SAVE
11
12  PUBLIC :: set_bounds
13
14CONTAINS
15
16  SUBROUTINE set_bounds_primal(d, cells, all, own)
17    TYPE(t_domain) :: d
18    TYPE(t_cellset) :: cells
19    LOGICAL :: all, own(:,:)         ! if all is .TRUE., include halo cells
20    REAL :: lon,lat
21    INTEGER :: i,j,k, n, halo_size
22
23    halo_size = MERGE(1,0,all)
24
25    ! count primal cells
26    n=0
27    DO j=d%jj_begin-halo_size, d%jj_end+halo_size
28       DO i=d%ii_begin-halo_size, d%ii_end+halo_size
29          IF (own(i,j) .OR. all ) n=n+1
30       ENDDO
31    ENDDO
32    cells%ncell = n
33   
34    IF ( .NOT. ALLOCATED(cells%ij) ) THEN 
35        ! now set bounds
36        ALLOCATE(cells%ij(n), cells%lon(n), cells%lat(n), cells%ind_glo(n))
37        ALLOCATE(cells%bnds_lon(0:5,n), cells%bnds_lat(0:5,n))
38    END IF
39
40    n=0
41    DO j=d%jj_begin-halo_size, d%jj_end+halo_size
42       DO i=d%ii_begin-halo_size, d%ii_end+halo_size
43          IF (own(i,j) .OR. all) THEN
44             n=n+1
45             CALL xyz2lonlat(d%xyz(:,i,j), lon, lat)
46             cells%lon(n)=lon*180./Pi
47             cells%lat(n)=lat*180./Pi
48             DO k=0,5
49                CALL xyz2lonlat(d%vertex(:,k,i,j), lon, lat)
50                cells%bnds_lon(k,n)=lon*180./Pi
51                cells%bnds_lat(k,n)=lat*180./Pi
52             END DO
53             cells%ij(n)=d%iim*(j-1)+i
54             cells%ind_glo(n) = d%assign_cell_glo(i,j)-1
55          END IF
56       END DO
57    END DO
58
59!    PRINT *, 'set_bounds_primal', all, halo_size, cells%ncell
60  END SUBROUTINE set_bounds_primal
61
62  SUBROUTINE set_bounds_dual(d, cells)
63    USE metric, ONLY : vup, vdown
64    TYPE(t_domain) :: d
65    TYPE(t_cellset) :: cells
66    REAL :: lonc, latc, lon(0:2), lat(0:2)
67    INTEGER :: i,j,k,n
68
69    ! count dual cells
70    n=0
71    DO j=d%jj_begin+1,d%jj_end
72       DO i=d%ii_begin,d%ii_end-1
73          n=n+2
74       ENDDO
75    ENDDO
76    cells%ncell = n
77
78    IF ( .NOT. ALLOCATED( cells%ind_glo ) ) THEN 
79        ! now set bounds
80        ALLOCATE(cells%ind_glo(n)) ! not set but must be allocated
81        ALLOCATE(cells%ij(n), cells%lon(n), cells%lat(n))
82        ALLOCATE(cells%bnds_lon(0:2,n), cells%bnds_lat(0:2,n))
83    END IF
84
85    n=0
86    DO j=d%jj_begin+1,d%jj_end
87       DO i=d%ii_begin,d%ii_end-1
88          n=n+1
89          CALL xyz2lonlat(d%vertex(:,vdown,i,j), lonc, latc)
90          CALL xyz2lonlat(d%xyz(:,i,j),     lon(0), lat(0))
91          CALL xyz2lonlat(d%xyz(:,i,j-1),   lon(1), lat(1))
92          CALL xyz2lonlat(d%xyz(:,i+1,j-1), lon(2), lat(2))
93          cells%lon(n)=lonc*180./Pi
94          cells%lat(n)=latc*180/Pi
95          DO k=0,2
96             cells%bnds_lat(k,n)=lat(k)*180./Pi
97             cells%bnds_lon(k,n)=lon(k)*180./Pi
98          END DO
99          cells%ij(n) = d%z_down + d%iim*(j-1)+i
100       ENDDO
101    ENDDO
102   
103    DO j=d%jj_begin,d%jj_end-1
104       DO i=d%ii_begin+1,d%ii_end
105          n=n+1
106          CALL xyz2lonlat(d%vertex(:,vup,i,j), lonc, latc)
107          CALL xyz2lonlat(d%xyz(:,i,j),     lon(0), lat(0))
108          CALL xyz2lonlat(d%xyz(:,i,j+1),   lon(1), lat(1))
109          CALL xyz2lonlat(d%xyz(:,i-1,j+1), lon(2), lat(2))
110          cells%lon(n)=lonc*180./Pi
111          cells%lat(n)=latc*180/Pi
112          DO k=0,2
113             cells%bnds_lat(k,n)=lat(k)*180./Pi
114             cells%bnds_lon(k,n)=lon(k)*180./Pi
115          END DO
116          cells%ij(n) = d%z_up + d%iim*(j-1)+i
117       ENDDO
118    ENDDO
119
120  END SUBROUTINE set_bounds_dual
121
122  SUBROUTINE set_bounds_edge(ind, d, cells)
123    USE metric, ONLY : cell_glo
124    INTEGER, INTENT(IN) :: ind
125    TYPE(t_domain) :: d
126    TYPE(t_cellset) :: cells
127    REAL :: lon(2), lat(2)
128    INTEGER :: i,j,ij,k,kk,n
129   
130    ! count edges
131    n=0
132    DO j=d%jj_begin,d%jj_end
133       DO i=d%ii_begin,d%ii_end
134          DO k=0,5
135             IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) &
136                  .AND. d%edge_assign_i(k,i,j)==i &
137                  .AND. d%edge_assign_j(k,i,j)==j &
138                  .AND. d%edge_assign_pos(k,i,j)==k) n=n+1
139          END DO
140       END DO
141    END DO
142    cells%ncell = n
143
144    IF ( .NOT. ALLOCATED( cells%ij ) ) THEN 
145        ! now set bounds
146        ALLOCATE(cells%ij(n), cells%lon(n), cells%lat(n), cells%ind_glo(n))
147        ALLOCATE(cells%sgn(n)) ! flip sign when reading/writing
148        ALLOCATE(cells%bnds_lon(2,n), cells%bnds_lat(2,n))
149    END IF
150
151    CALL swap_dimensions(ind)
152    CALL swap_geometry(ind)
153   
154    n=0
155    DO j=d%jj_begin,d%jj_end
156       DO i=d%ii_begin,d%ii_end
157          DO k=0,5
158             IF (d%edge_assign_domain(k,i,j)==domloc_glo_ind(ind) &
159                  .AND. d%edge_assign_i(k,i,j)==i &
160                  .AND. d%edge_assign_j(k,i,j)==j &
161                  .AND. d%edge_assign_pos(k,i,j)==k) THEN
162                n=n+1
163                ij=(j-1)*d%iim+i+u_pos(k+1)
164                kk = MOD(k+d%delta(i,j)+6,6)
165                cells%ij(n)     = ij
166                cells%sgn(n)    = d%edge_assign_sign(k,i,j)
167                cells%ind_glo(n)= cell_glo(d%assign_cell_glo(i,j))%edge(kk)-1
168                cells%lon(n)    = lon_e(ij)*180./Pi
169                cells%lat(n)    = lat_e(ij)*180./Pi
170               
171                kk = MOD(k-1+6,6)
172                CALL xyz2lonlat(d%vertex(:,kk,i,j), lon(1),lat(1))
173                CALL xyz2lonlat(d%vertex(:,k, i,j), lon(2),lat(2))
174                cells%bnds_lon(:,n)=lon(:)*180./Pi
175                cells%bnds_lat(:,n)=lat(:)*180/Pi
176             END IF
177          END DO
178       END DO
179    END DO
180
181  END SUBROUTINE set_bounds_edge
182
183  SUBROUTINE set_bounds(domain_type, glo)
184    TYPE(t_domain), POINTER :: domain_type(:), d
185    LOGICAL :: glo
186    INTEGER :: ind
187    !$OMP BARRIER   
188    !$OMP MASTER
189    ! IF glo is .TRUE. we are dealing with the global mesh, otherwise with the local mesh
190    ! write_field uses the global mesh and may want halo cells
191    ! output_field (XIOS) uses the local mesh and uses only own cells
192    ! output_field uses edges while write_field uses only primal and dual cells
193    DO ind=1, SIZE(domain_type)
194       d=>domain_type(ind)
195       IF(glo) THEN ! global mesh / write_field
196          ! primal cell i,j is owned if d%assign_domain(i,j)==ind
197          CALL set_bounds_primal(d, d%primal_own, .FALSE., d%assign_domain==ind)
198          CALL set_bounds_primal(d, d%primal_all, .TRUE., d%assign_domain==ind)
199          CALL set_bounds_dual(d, d%dual_own)
200          CALL set_bounds_dual(d, d%dual_all)
201       ELSE  ! local mesh / XIOS
202          ! primal cell i,j is owned if d%own(i,j)==.TRUE.
203          CALL set_bounds_primal(d, d%primal_own, .FALSE., d%own)
204          CALL set_bounds_dual(d, d%dual_own)
205          CALL set_bounds_edge(ind, d, d%edge_own)
206       END IF
207    END DO
208    !$OMP END MASTER
209    !$OMP BARRIER
210  END SUBROUTINE set_bounds
211 
212END MODULE set_bounds_mod
Note: See TracBrowser for help on using the repository browser.