source: codes/icosagcm/devel/src/sphere/geometry.f90 @ 879

Last change on this file since 879 was 879, checked in by dubos, 5 years ago

devel : introduced derived type to store cell bounds

File size: 6.7 KB
Line 
1MODULE geometry
2  USE field_mod
3  IMPLICIT NONE
4
5  TYPE t_geometry
6    TYPE(t_field),POINTER :: centroid(:)
7    TYPE(t_field),POINTER :: xyz_i(:)
8    TYPE(t_field),POINTER :: xyz_e(:)
9    TYPE(t_field),POINTER :: xyz_v(:)
10    TYPE(t_field),POINTER :: lon_i(:)
11    TYPE(t_field),POINTER :: lon_e(:)
12    TYPE(t_field),POINTER :: lat_i(:)
13    TYPE(t_field),POINTER :: lat_e(:)
14    TYPE(t_field),POINTER :: ep_e(:)
15    TYPE(t_field),POINTER :: et_e(:)
16    TYPE(t_field),POINTER :: elon_i(:)
17    TYPE(t_field),POINTER :: elat_i(:)
18    TYPE(t_field),POINTER :: elon_e(:)
19    TYPE(t_field),POINTER :: elat_e(:)
20    TYPE(t_field),POINTER :: Ai(:)
21    TYPE(t_field),POINTER :: Av(:)
22    TYPE(t_field),POINTER :: de(:)
23    TYPE(t_field),POINTER :: le(:)
24    TYPE(t_field),POINTER :: le_de(:) ! le/de, 0. if de=0.
25    TYPE(t_field),POINTER :: Riv(:)
26    TYPE(t_field),POINTER :: S1(:)
27    TYPE(t_field),POINTER :: S2(:)
28    TYPE(t_field),POINTER :: Riv2(:)
29    TYPE(t_field),POINTER :: ne(:)
30    TYPE(t_field),POINTER :: Wee(:)
31    TYPE(t_field),POINTER :: bi(:)
32    TYPE(t_field),POINTER :: fv(:)
33  END TYPE t_geometry   
34
35  TYPE(t_geometry),SAVE,TARGET :: geom
36
37 
38  REAL(rstd),POINTER :: Ai(:)          ! area of a cell
39!$OMP THREADPRIVATE(Ai)
40  REAL(rstd),POINTER :: centroid(:,:)  ! coordinate of the centroid of the cell
41!$OMP THREADPRIVATE(centroid)
42  REAL(rstd),POINTER :: xyz_i(:,:)     ! coordinate of the center of the cell (voronoi)
43!$OMP THREADPRIVATE(xyz_i)
44  REAL(rstd),POINTER :: xyz_e(:,:)     ! coordinate of a wind point on the cell on a edge
45!$OMP THREADPRIVATE(xyz_e)
46  REAL(rstd),POINTER :: xyz_v(:,:)     ! coordinate of a vertex (center of the dual mesh)
47!$OMP THREADPRIVATE(xyz_v)
48  REAL(rstd),POINTER :: lon_i(:)       ! longitude of the center of the cell (voronoi)
49!$OMP THREADPRIVATE(lon_i)
50  REAL(rstd),POINTER :: lon_e(:)       ! longitude of a wind point on the cell on a edge
51!$OMP THREADPRIVATE(lon_e)
52  REAL(rstd),POINTER :: lat_i(:)       ! latitude of the center of the cell (voronoi)
53!$OMP THREADPRIVATE(lat_i)
54  REAL(rstd),POINTER :: lat_e(:)       ! latitude of a wind point on the cell on a edge
55!$OMP THREADPRIVATE(lat_e)
56  REAL(rstd),POINTER :: ep_e(:,:)      ! perpendicular unit vector of a edge (outsider)
57!$OMP THREADPRIVATE(ep_e)
58  REAL(rstd),POINTER :: et_e(:,:)      ! tangeantial unit vector of a edge
59!$OMP THREADPRIVATE(et_e)
60  REAL(rstd),POINTER :: elon_i(:,:)    ! unit longitude vector on the center
61!$OMP THREADPRIVATE(elon_i)
62  REAL(rstd),POINTER :: elat_i(:,:)    ! unit latitude vector on the center
63!$OMP THREADPRIVATE(elat_i)
64  REAL(rstd),POINTER :: elon_e(:,:)    ! unit longitude vector on a wind point
65!$OMP THREADPRIVATE(elon_e)
66  REAL(rstd),POINTER :: elat_e(:,:)    ! unit latitude vector on a wind point
67!$OMP THREADPRIVATE(elat_e)
68  REAL(rstd),POINTER :: Av(:)          ! area of dual mesk cell
69!$OMP THREADPRIVATE(Av)
70  REAL(rstd),POINTER :: de(:)          ! distance from a neighbour == lenght of an edge of the dual mesh
71!$OMP THREADPRIVATE(de)
72  REAL(rstd),POINTER :: le(:)          ! lenght of a edge
73!$OMP THREADPRIVATE(le)
74  REAL(rstd),POINTER :: le_de(:)       ! le/de
75!$OMP THREADPRIVATE(le_de)
76  REAL(rstd),POINTER :: S1(:,:)        ! area of sub-triangle
77!$OMP THREADPRIVATE(S1)
78  REAL(rstd),POINTER :: S2(:,:)        ! area of sub-tirangle
79!$OMP THREADPRIVATE(S2)
80  REAL(rstd),POINTER :: Riv(:,:)       ! weight
81!$OMP THREADPRIVATE(Riv)
82  REAL(rstd),POINTER :: Riv2(:,:)      ! weight
83!$OMP THREADPRIVATE(Riv2)
84  INTEGER,POINTER    :: ne(:,:)        ! convention for the way on the normal wind on an edge
85!$OMP THREADPRIVATE(ne)
86  REAL(rstd),POINTER :: Wee(:,:,:)     ! weight
87!$OMP THREADPRIVATE(Wee)
88  REAL(rstd),POINTER :: bi(:)          ! orographie
89!$OMP THREADPRIVATE(bi)
90  REAL(rstd),POINTER :: fv(:)          ! coriolis (evaluted on a vertex)
91!$OMP THREADPRIVATE(fv)
92     
93
94  INTEGER, PARAMETER :: ne_right=1
95  INTEGER, PARAMETER :: ne_rup=-1
96  INTEGER, PARAMETER :: ne_lup=1
97  INTEGER, PARAMETER :: ne_left=-1
98  INTEGER, PARAMETER :: ne_ldown=1
99  INTEGER, PARAMETER :: ne_rdown=-1
100
101CONTAINS
102
103  SUBROUTINE allocate_geometry
104  USE field_mod
105  IMPLICIT NONE
106  INTEGER, PARAMETER :: nvertex=6    ! FIXME unstructured
107
108    CALL allocate_field(geom%Ai,field_t,type_real,name='Ai')
109
110    CALL allocate_field(geom%xyz_i,field_t,type_real,3)
111    CALL allocate_field(geom%lon_i,field_t,type_real, name='lon_i')
112    CALL allocate_field(geom%lat_i,field_t,type_real, name='lat_i')
113    CALL allocate_field(geom%elon_i,field_t,type_real,3)
114    CALL allocate_field(geom%elat_i,field_t,type_real,3)
115    CALL allocate_field(geom%centroid,field_t,type_real,3)
116
117    CALL allocate_field(geom%xyz_e,field_u,type_real,3)
118    CALL allocate_field(geom%lon_e,field_u,type_real, name='lon_e')
119    CALL allocate_field(geom%lat_e,field_u,type_real, name='lat_e')
120    CALL allocate_field(geom%elon_e,field_u,type_real,3)
121    CALL allocate_field(geom%elat_e,field_u,type_real,3)
122    CALL allocate_field(geom%ep_e,field_u,type_real,3)
123    CALL allocate_field(geom%et_e,field_u,type_real,3)
124
125    CALL allocate_field(geom%xyz_v,field_z,type_real,3)
126    CALL allocate_field(geom%de,field_u,type_real, name='de')
127    CALL allocate_field(geom%le,field_u,type_real, name='le')
128    CALL allocate_field(geom%le_de,field_u,type_real, name='le_de')
129    CALL allocate_field(geom%bi,field_t,type_real)
130    CALL allocate_field(geom%Av,field_z,type_real, name='Av')
131    CALL allocate_field(geom%S1,field_t,type_real,nvertex)
132    CALL allocate_field(geom%S2,field_t,type_real,nvertex)
133    CALL allocate_field(geom%Riv,field_t,type_real,nvertex)
134    CALL allocate_field(geom%Riv2,field_t,type_real,nvertex)
135    CALL allocate_field(geom%ne,field_t,type_integer,nvertex)
136    CALL allocate_field(geom%Wee,field_u,type_real,5,2) ! FIXME unstructured
137    CALL allocate_field(geom%bi,field_t,type_real)
138    CALL allocate_field(geom%fv,field_z,type_real, name='fv')
139
140  END SUBROUTINE allocate_geometry
141 
142 
143  SUBROUTINE swap_geometry(ind)
144  USE field_mod
145  USE domain_mod, ONLY : swap_needed
146  IMPLICIT NONE
147    INTEGER,INTENT(IN) :: ind
148
149    IF(.NOT. swap_needed) RETURN
150
151!!$OMP MASTER
152    Ai=geom%Ai(ind)
153    xyz_i=geom%xyz_i(ind)
154    centroid=geom%centroid(ind)
155    xyz_e=geom%xyz_e(ind)
156    ep_e=geom%ep_e(ind)
157    et_e=geom%et_e(ind)
158    lon_i=geom%lon_i(ind)
159    lat_i=geom%lat_i(ind)
160    lon_e=geom%lon_e(ind)
161    lat_e=geom%lat_e(ind)
162    elon_i=geom%elon_i(ind)
163    elat_i=geom%elat_i(ind)
164    elon_e=geom%elon_e(ind)
165    elat_e=geom%elat_e(ind)
166    xyz_v=geom%xyz_v(ind)
167    de=geom%de(ind)
168    le=geom%le(ind)
169    le_de=geom%le_de(ind)
170    Av=geom%Av(ind)
171    S1=geom%S1(ind)
172    S2=geom%S2(ind)
173    Riv=geom%Riv(ind)
174    Riv2=geom%Riv2(ind)
175    ne=geom%ne(ind)
176    Wee=geom%Wee(ind)
177    bi=geom%bi(ind)
178    fv=geom%fv(ind)
179!!$OMP END MASTER
180!!$OMP BARRIER   
181  END SUBROUTINE swap_geometry
182
183END MODULE geometry
Note: See TracBrowser for help on using the repository browser.