source: codes/icosagcm/trunk/src/geometry.f90 @ 356

Last change on this file since 356 was 356, checked in by dubos, 9 years ago

Moved Schmidt transform after grid optimization ; pass correct mesh info to XIOS
File size: 20.5 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 :: Riv(:)
25    TYPE(t_field),POINTER :: Riv2(:)
26    TYPE(t_field),POINTER :: ne(:)
27    TYPE(t_field),POINTER :: Wee(:)
28    TYPE(t_field),POINTER :: bi(:)
29    TYPE(t_field),POINTER :: fv(:)
30  END TYPE t_geometry   
31
32  TYPE(t_geometry),SAVE,TARGET :: geom
33
34 
35  REAL(rstd),POINTER :: Ai(:)          ! area of a cell
36!$OMP THREADPRIVATE(Ai)
37  REAL(rstd),POINTER :: centroid(:,:)  ! coordinate of the centroid of the cell
38!$OMP THREADPRIVATE(centroid)
39  REAL(rstd),POINTER :: xyz_i(:,:)     ! coordinate of the center of the cell (voronoi)
40!$OMP THREADPRIVATE(xyz_i)
41  REAL(rstd),POINTER :: xyz_e(:,:)     ! coordinate of a wind point on the cell on a edge
42!$OMP THREADPRIVATE(xyz_e)
43  REAL(rstd),POINTER :: xyz_v(:,:)     ! coordinate of a vertex (center of the dual mesh)
44!$OMP THREADPRIVATE(xyz_v)
45  REAL(rstd),POINTER :: lon_i(:)       ! longitude of the center of the cell (voronoi)
46!$OMP THREADPRIVATE(lon_i)
47  REAL(rstd),POINTER :: lon_e(:)       ! longitude of a wind point on the cell on a edge
48!$OMP THREADPRIVATE(lon_e)
49  REAL(rstd),POINTER :: lat_i(:)       ! latitude of the center of the cell (voronoi)
50!$OMP THREADPRIVATE(lat_i)
51  REAL(rstd),POINTER :: lat_e(:)       ! latitude of a wind point on the cell on a edge
52!$OMP THREADPRIVATE(lat_e)
53  REAL(rstd),POINTER :: ep_e(:,:)      ! perpendicular unit vector of a edge (outsider)
54!$OMP THREADPRIVATE(ep_e)
55  REAL(rstd),POINTER :: et_e(:,:)      ! tangeantial unit vector of a edge
56!$OMP THREADPRIVATE(et_e)
57  REAL(rstd),POINTER :: elon_i(:,:)    ! unit longitude vector on the center
58!$OMP THREADPRIVATE(elon_i)
59  REAL(rstd),POINTER :: elat_i(:,:)    ! unit latitude vector on the center
60!$OMP THREADPRIVATE(elat_i)
61  REAL(rstd),POINTER :: elon_e(:,:)    ! unit longitude vector on a wind point
62!$OMP THREADPRIVATE(elon_e)
63  REAL(rstd),POINTER :: elat_e(:,:)    ! unit latitude vector on a wind point
64!$OMP THREADPRIVATE(elat_e)
65  REAL(rstd),POINTER :: Av(:)          ! area of dual mesk cell
66!$OMP THREADPRIVATE(Av)
67  REAL(rstd),POINTER :: de(:)          ! distance from a neighbour == lenght of an edge of the dual mesh
68!$OMP THREADPRIVATE(de)
69  REAL(rstd),POINTER :: le(:)          ! lenght of a edge
70!$OMP THREADPRIVATE(le)
71  REAL(rstd),POINTER :: Riv(:,:)       ! weight
72!$OMP THREADPRIVATE(Riv)
73  REAL(rstd),POINTER :: Riv2(:,:)      ! weight
74!$OMP THREADPRIVATE(Riv2)
75  INTEGER,POINTER    :: ne(:,:)        ! convention for the way on the normal wind on an edge
76!$OMP THREADPRIVATE(ne)
77  REAL(rstd),POINTER :: Wee(:,:,:)     ! weight
78!$OMP THREADPRIVATE(Wee)
79  REAL(rstd),POINTER :: bi(:)          ! orographie
80!$OMP THREADPRIVATE(bi)
81  REAL(rstd),POINTER :: fv(:)          ! coriolis (evaluted on a vertex)
82!$OMP THREADPRIVATE(fv)
83     
84
85  INTEGER, PARAMETER :: ne_right=1
86  INTEGER, PARAMETER :: ne_rup=-1
87  INTEGER, PARAMETER :: ne_lup=1
88  INTEGER, PARAMETER :: ne_left=-1
89  INTEGER, PARAMETER :: ne_ldown=1
90  INTEGER, PARAMETER :: ne_rdown=-1
91
92CONTAINS
93
94  SUBROUTINE allocate_geometry
95  USE field_mod
96  IMPLICIT NONE
97 
98    CALL allocate_field(geom%Ai,field_t,type_real,name='Ai')
99
100    CALL allocate_field(geom%xyz_i,field_t,type_real,3)
101    CALL allocate_field(geom%lon_i,field_t,type_real)
102    CALL allocate_field(geom%lat_i,field_t,type_real)
103    CALL allocate_field(geom%elon_i,field_t,type_real,3)
104    CALL allocate_field(geom%elat_i,field_t,type_real,3)
105    CALL allocate_field(geom%centroid,field_t,type_real,3)
106
107    CALL allocate_field(geom%xyz_e,field_u,type_real,3)
108    CALL allocate_field(geom%lon_e,field_u,type_real)
109    CALL allocate_field(geom%lat_e,field_u,type_real)
110    CALL allocate_field(geom%elon_e,field_u,type_real,3)
111    CALL allocate_field(geom%elat_e,field_u,type_real,3)
112    CALL allocate_field(geom%ep_e,field_u,type_real,3)
113    CALL allocate_field(geom%et_e,field_u,type_real,3)
114
115    CALL allocate_field(geom%xyz_v,field_z,type_real,3)
116    CALL allocate_field(geom%de,field_u,type_real)
117    CALL allocate_field(geom%le,field_u,type_real)
118    CALL allocate_field(geom%bi,field_t,type_real)
119    CALL allocate_field(geom%Av,field_z,type_real)
120    CALL allocate_field(geom%Riv,field_t,type_real,6)
121    CALL allocate_field(geom%Riv2,field_t,type_real,6)
122    CALL allocate_field(geom%ne,field_t,type_integer,6)
123    CALL allocate_field(geom%Wee,field_u,type_real,5,2)
124    CALL allocate_field(geom%bi,field_t,type_real)
125    CALL allocate_field(geom%fv,field_z,type_real)
126
127  END SUBROUTINE allocate_geometry
128 
129 
130  SUBROUTINE swap_geometry(ind)
131  USE field_mod
132  IMPLICIT NONE
133    INTEGER,INTENT(IN) :: ind
134!!$OMP MASTER
135    Ai=geom%Ai(ind)
136    xyz_i=geom%xyz_i(ind)
137    centroid=geom%centroid(ind)
138    xyz_e=geom%xyz_e(ind)
139    ep_e=geom%ep_e(ind)
140    et_e=geom%et_e(ind)
141    lon_i=geom%lon_i(ind)
142    lat_i=geom%lat_i(ind)
143    lon_e=geom%lon_e(ind)
144    lat_e=geom%lat_e(ind)
145    elon_i=geom%elon_i(ind)
146    elat_i=geom%elat_i(ind)
147    elon_e=geom%elon_e(ind)
148    elat_e=geom%elat_e(ind)
149    xyz_v=geom%xyz_v(ind)
150    de=geom%de(ind)
151    le=geom%le(ind)
152    Av=geom%Av(ind)
153    Riv=geom%Riv(ind)
154    Riv2=geom%Riv2(ind)
155    ne=geom%ne(ind)
156    Wee=geom%Wee(ind)
157    bi=geom%bi(ind)
158    fv=geom%fv(ind)
159!!$OMP END MASTER
160!!$OMP BARRIER   
161  END SUBROUTINE swap_geometry
162
163  SUBROUTINE update_circumcenters 
164    USE domain_mod
165    USE dimensions
166    USE spherical_geom_mod
167    USE vector
168    USE transfert_mod
169    USE omp_para
170
171    IMPLICIT NONE
172    REAL(rstd) :: x1(3),x2(3)
173    REAL(rstd) :: vect(3,6)
174    REAL(rstd) :: centr(3)
175    INTEGER :: ind,i,j,n,k
176    TYPE(t_message),SAVE :: message0, message1
177    LOGICAL, SAVE :: first=.TRUE.
178!$OMP THREADPRIVATE(first)   
179   
180    IF (first) THEN
181      CALL init_message(geom%xyz_i, req_i0 ,message0)
182      CALL init_message(geom%xyz_i, req_i1 ,message1)
183      first=.FALSE.
184    ENDIF
185   
186    CALL transfert_message(geom%xyz_i,message0)
187    CALL transfert_message(geom%xyz_i,message1)
188   
189    DO ind=1,ndomain
190      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
191      CALL swap_dimensions(ind)
192      CALL swap_geometry(ind)
193      DO j=jj_begin,jj_end
194        DO i=ii_begin,ii_end
195          n=(j-1)*iim+i
196          DO k=0,5
197            x1(:) = xyz_i(n+t_pos(k+1),:)
198            x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:)
199            if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:)
200            CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:))
201          ENDDO
202        ENDDO
203      ENDDO
204    ENDDO
205
206  END SUBROUTINE update_circumcenters
207
208  SUBROUTINE remap_schmidt_loc
209    USE spherical_geom_mod
210    USE getin_mod
211    USE omp_para
212    USE domain_mod
213    USE dimensions
214    IMPLICIT NONE
215    INTEGER :: ind,i,j,n
216    REAL(rstd) :: schmidt_factor, schmidt_lon, schmidt_lat
217
218    ! Schmidt transform parameters
219    schmidt_factor = 1.
220    CALL getin('schmidt_factor', schmidt_factor)
221    schmidt_factor =  schmidt_factor**2.
222   
223    schmidt_lon = 0.
224    CALL getin('schmidt_lon', schmidt_lon)
225    schmidt_lon = schmidt_lon * pi/180.
226
227    schmidt_lat = 45.
228    CALL getin('schmidt_lat', schmidt_lat)
229    schmidt_lat = schmidt_lat * pi/180.
230
231    DO ind=1,ndomain
232      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
233      CALL swap_dimensions(ind)
234      CALL swap_geometry(ind)
235      DO j=jj_begin,jj_end
236        DO i=ii_begin,ii_end
237          n=(j-1)*iim+i
238          CALL schmidt_transform(xyz_i(n,:), schmidt_factor, schmidt_lon, schmidt_lat)
239        ENDDO
240      ENDDO
241    ENDDO
242  END SUBROUTINE remap_schmidt_loc
243
244  SUBROUTINE optimize_geometry
245  USE metric
246  USE spherical_geom_mod
247  USE domain_mod
248  USE dimensions
249  USE transfert_mod
250  USE vector
251  USE getin_mod
252  USE omp_para
253  IMPLICIT NONE
254    INTEGER :: nb_it=0
255    TYPE(t_domain),POINTER :: d
256    INTEGER :: ind,it,i,j,n,k
257    REAL(rstd) :: x1(3),x2(3)
258    REAL(rstd) :: vect(3,6)
259    REAL(rstd) :: centr(3)
260    REAL(rstd) :: sum
261    LOGICAL    :: check
262   
263   
264    CALL getin('optim_it',nb_it)
265   
266    DO ind=1,ndomain
267      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master) CYCLE
268      d=>domain(ind)
269      CALL swap_dimensions(ind)
270      CALL swap_geometry(ind)
271      DO j=jj_begin,jj_end
272        DO i=ii_begin,ii_end
273          n=(j-1)*iim+i
274          xyz_i(n,:)=d%xyz(:,i,j) 
275        ENDDO
276      ENDDO
277    ENDDO
278   
279    CALL update_circumcenters
280
281    DO ind=1,ndomain
282      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
283      d=>domain(ind)
284      CALL swap_dimensions(ind)
285      CALL swap_geometry(ind)
286      DO j=jj_begin,jj_end
287        DO i=ii_begin,ii_end
288          n=(j-1)*iim+i
289          DO k=0,5
290            x1(:) = xyz_v(n+z_pos(k+1),:)
291            x2(:) = d%vertex(:,k,i,j) 
292            IF (norm(x1-x2)>1e-10) THEN
293              PRINT*,"vertex diff ",ind,i,j,k
294              PRINT*,x1
295              PRINT*,x2
296            ENDIF
297          ENDDO
298        ENDDO
299      ENDDO
300    ENDDO
301   
302   
303    DO it=1,nb_it
304      IF (MOD(it,100)==0) THEN
305        check=is_master
306      ELSE
307        check=.FALSE.
308      ENDIF
309     
310      sum=0
311      DO ind=1,ndomain
312      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
313        CALL swap_dimensions(ind)
314        CALL swap_geometry(ind)
315        DO j=jj_begin,jj_end
316          DO i=ii_begin,ii_end
317            n=(j-1)*iim+i
318            vect(:,1)=xyz_v(n+z_rup,:)
319            vect(:,2)=xyz_v(n+z_up,:)
320            vect(:,3)=xyz_v(n+z_lup,:)
321            vect(:,4)=xyz_v(n+z_ldown,:)
322            vect(:,5)=xyz_v(n+z_down,:)
323            vect(:,6)=xyz_v(n+z_rdown,:)
324            CALL compute_centroid(vect,6,centr)
325            IF (check) THEN
326              sum=MAX(sum,norm(xyz_i(n,:)-centr(:)))
327            ENDIF
328            xyz_i(n,:)=centr(:)
329          ENDDO
330        ENDDO
331       
332      ENDDO
333
334      IF (check) PRINT *,"it = ",it,"  diff centroid circumcenter ",sum
335
336     CALL update_circumcenters
337
338    ENDDO
339   
340  END SUBROUTINE optimize_geometry
341
342  SUBROUTINE update_domain 
343    ! copy position of generators and vertices back into domain(:)%xyz/vertex
344    ! so that XIOS gets the right values
345    USE omp_para
346    USE dimensions
347    USE domain_mod
348    INTEGER :: ind,i,j,k,n
349    TYPE(t_domain),POINTER :: d
350    DO ind=1,ndomain
351       IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
352       d=>domain(ind)
353       CALL swap_dimensions(ind)
354       CALL swap_geometry(ind)
355       DO j=jj_begin,jj_end
356          DO i=ii_begin,ii_end         
357             n=(j-1)*iim+i
358             d%xyz(:,i,j)=xyz_i(n,:)
359             DO k=0,5
360                d%vertex(:,k,i,j) = xyz_v(n+z_pos(k+1),:)
361             END DO
362          END DO
363       END DO
364    END DO
365  END SUBROUTINE update_domain
366 
367  SUBROUTINE set_geometry
368  USE metric
369  USE vector
370  USE spherical_geom_mod
371  USE domain_mod
372  USE dimensions
373  USE transfert_mod
374  USE getin_mod
375  USE omp_para
376  IMPLICIT NONE
377
378    REAL(rstd) :: surf(6)
379    REAL(rstd) :: surf_v(6)
380    REAL(rstd) :: vect(3,6)
381    REAL(rstd) :: centr(3)
382    REAL(rstd) :: vet(3),vep(3), vertex(3)
383    INTEGER :: ind,i,j,k,n
384    TYPE(t_domain),POINTER :: d
385    REAL(rstd) ::  S
386    REAL(rstd) :: w(6)
387    REAL(rstd) :: lon,lat
388    INTEGER :: ii_glo,jj_glo
389    REAL(rstd) :: S1,S2
390         
391     
392    CALL optimize_geometry
393    CALL remap_schmidt_loc
394    CALL update_circumcenters
395    ! copy position of generators and vertices back into domain(:)%xyz/vertex
396    ! so that XIOS gets the right values
397    CALL update_domain
398
399    DO ind=1,ndomain
400      IF (.NOT. assigned_domain(ind)  .OR. .NOT. is_omp_level_master ) CYCLE
401      d=>domain(ind)
402      CALL swap_dimensions(ind)
403      CALL swap_geometry(ind)
404      DO j=jj_begin-1,jj_end+1
405        DO i=ii_begin-1,ii_end+1
406          n=(j-1)*iim+i
407
408          DO k=0,5
409            ne(n,k+1)=d%ne(k,i,j)
410          ENDDO
411
412          vect(:,1)=xyz_v(n+z_rup,:)
413          vect(:,2)=xyz_v(n+z_up,:)
414          vect(:,3)=xyz_v(n+z_lup,:)
415          vect(:,4)=xyz_v(n+z_ldown,:)
416          vect(:,5)=xyz_v(n+z_down,:)
417          vect(:,6)=xyz_v(n+z_rdown,:)
418          CALL compute_centroid(vect,6,centr)
419          centroid(n,:)=centr(:)
420
421               
422          CALL xyz2lonlat(xyz_v(n+z_up,:),lon,lat)
423          fv(n+z_up)=2*sin(lat)*omega
424          CALL xyz2lonlat(xyz_v(n+z_down,:),lon,lat)
425          fv(n+z_down)=2*sin(lat)*omega
426         
427          bi(n)=0. 
428       
429          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_right,:),de(n+u_right))
430          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_lup,:),de(n+u_lup))
431          CALL dist_cart(xyz_i(n,:),xyz_i(n+t_ldown,:),de(n+u_ldown))
432         
433          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_right,:),0.5,xyz_e(n+u_right,:))
434          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_lup,:),0.5,xyz_e(n+u_lup,:))
435          CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_ldown,:),0.5,xyz_e(n+u_ldown,:))
436
437          CALL dist_cart(xyz_v(n+z_rdown,:), xyz_v(n+z_rup,:),le(n+u_right))
438          CALL dist_cart(xyz_v(n+z_up,:), xyz_v(n+z_lup,:),le(n+u_lup))
439          CALL dist_cart(xyz_v(n+z_ldown,:), xyz_v(n+z_down,:),le(n+u_ldown))
440         
441          Ai(n)=0
442          DO k=0,5
443            CALL surf_triangle(xyz_i(n,:),xyz_i(n+t_pos(k+1),:),xyz_i(n+t_pos(MOD((k+1+6),6)+1),:),surf_v(k+1))
444            CALL surf_triangle(xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:),surf(k+1))
445            Ai(n)=Ai(n)+surf(k+1)
446            IF (i==ii_end .AND. j==jj_begin) THEN
447              IF (Ai(n)<1e20) THEN
448              ELSE
449                PRINT *,"PB !!",Ai(n),k,surf(k+1)
450                PRINT*,xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:)
451              ENDIF
452            ENDIF
453          ENDDO
454
455          ! Sign convention : Ringler et al., JCP 2010, eq. 21 p. 3071
456          ! Normal component is along outgoing normal vector if ne=1
457
458          CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),vep)
459          IF (norm(vep)>1e-30) THEN
460            vep(:)=vep(:)/norm(vep)                         ! Inward normal vector
461            CALL cross_product2(vep,xyz_e(n+u_right,:),vet) ! Counter-clockwise tangent vector
462            vet(:)=vet(:)/norm(vet)
463            ep_e(n+u_right,:)=-vep(:)*ne(n,right)
464            et_e(n+u_right,:)=vet(:)*ne(n,right)
465          ENDIF
466
467          CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),vep)
468          IF (norm(vep)>1e-30) THEN
469            vep(:)=vep(:)/norm(vep)
470            CALL cross_product2(vep,xyz_e(n+u_lup,:),vet)
471            vet(:)=vet(:)/norm(vet)
472            ep_e(n+u_lup,:)=-vep(:)*ne(n,lup)
473            et_e(n+u_lup,:)=vet(:)*ne(n,lup)
474          ENDIF
475
476          CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),vep)
477          IF (norm(vep)>1e-30) THEN
478            vep(:)=vep(:)/norm(vep)
479            CALL cross_product2(vep,xyz_e(n+u_ldown,:),vet)
480            vet(:)=vet(:)/norm(vet)
481            ep_e(n+u_ldown,:)=-vep(:)*ne(n,ldown)
482            et_e(n+u_ldown,:)=vet(:)*ne(n,ldown)
483          ENDIF
484
485          CALL xyz2lonlat(xyz_i(n,:),lon,lat)
486          lon_i(n)=lon
487          lat_i(n)=lat
488          elon_i(n,1) = -sin(lon)
489          elon_i(n,2) = cos(lon)
490          elon_i(n,3) = 0
491          elat_i(n,1) = -cos(lon)*sin(lat)
492          elat_i(n,2) = -sin(lon)*sin(lat)
493          elat_i(n,3) = cos(lat)
494
495         
496          CALL xyz2lonlat(xyz_e(n+u_right,:),lon,lat)
497          lon_e(n+u_right)=lon
498          lat_e(n+u_right)=lat
499          elon_e(n+u_right,1) = -sin(lon)
500          elon_e(n+u_right,2) = cos(lon)
501          elon_e(n+u_right,3) = 0
502          elat_e(n+u_right,1) = -cos(lon)*sin(lat)
503          elat_e(n+u_right,2) = -sin(lon)*sin(lat)
504          elat_e(n+u_right,3) = cos(lat)
505         
506          CALL xyz2lonlat(xyz_e(n+u_lup,:),lon,lat)
507          lon_e(n+u_lup)=lon
508          lat_e(n+u_lup)=lat
509          elon_e(n+u_lup,1) = -sin(lon)
510          elon_e(n+u_lup,2) = cos(lon)
511          elon_e(n+u_lup,3) = 0
512          elat_e(n+u_lup,1) = -cos(lon)*sin(lat)
513          elat_e(n+u_lup,2) = -sin(lon)*sin(lat)
514          elat_e(n+u_lup,3) = cos(lat)
515 
516          CALL xyz2lonlat(xyz_e(n+u_ldown,:),lon,lat)
517          lon_e(n+u_ldown)=lon
518          lat_e(n+u_ldown)=lat
519          elon_e(n+u_ldown,1) = -sin(lon)
520          elon_e(n+u_ldown,2) = cos(lon)
521          elon_e(n+u_ldown,3) = 0
522          elat_e(n+u_ldown,1) = -cos(lon)*sin(lat)
523          elat_e(n+u_ldown,2) = -sin(lon)*sin(lat)
524          elat_e(n+u_ldown,3) = cos(lat)
525 
526         
527          DO k=0,5
528            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1)
529            CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(MOD(k+1+6,6)+1),:),S2)
530            Riv(n,k+1)=0.5*(S1+S2)/Ai(n)
531            Riv2(n,k+1)=0.5*(S1+S2)/surf_v(k+1)
532          ENDDO
533
534          DO k=1,6
535            IF (ABS(surf_v(k))<1e-30) THEN
536              Riv(n,k)=0.
537            ENDIF
538          ENDDO
539         
540          Av(n+z_up)=surf_v(vup)+1e-100
541          Av(n+z_down)=surf_v(vdown)+1e-100
542       
543        ENDDO
544      ENDDO
545     
546      DO j=jj_begin,jj_end
547        DO i=ii_begin,ii_end
548          n=(j-1)*iim+i
549   
550          CALL compute_wee(n,right,w)
551          Wee(n+u_right,:,1)=w(1:5)
552
553          CALL compute_wee(n+t_right,left,w)
554          Wee(n+u_right,:,2)=w(1:5)
555
556
557          CALL compute_wee(n,lup,w)
558          Wee(n+u_lup,:,1)=w(1:5)
559
560          CALL compute_wee(n+t_lup,rdown,w)
561          Wee(n+u_lup,:,2)=w(1:5)
562
563
564          CALL compute_wee(n,ldown,w)
565          Wee(n+u_ldown,:,1)=w(1:5)
566
567          CALL compute_wee(n+t_ldown,rup,w)
568          Wee(n+u_ldown,:,2)=w(1:5)
569
570        ENDDO
571      ENDDO
572     
573      DO j=jj_begin,jj_end
574        DO i=ii_begin,ii_end
575          n=(j-1)*iim+i
576          ii_glo=d%ii_begin_glo-d%ii_begin+i
577          jj_glo=d%jj_begin_glo-d%jj_begin+j
578         
579          IF (ii_glo==1 .AND. jj_glo==1) THEN
580            le(n+u_ldown)=0
581            xyz_v(n+z_ldown,:)=xyz_v(n+z_down,:)
582                       
583          ENDIF
584
585          IF (ii_glo==iim_glo .AND. jj_glo==1) THEN
586            le(n+u_right)=0
587            xyz_v(n+z_rdown,:)=xyz_v(n+z_rup,:)
588          ENDIF
589
590          IF (ii_glo==iim_glo .AND. jj_glo==jjm_glo) THEN
591            le(n+u_rup)=0
592            xyz_v(n+z_rup,:)=xyz_v(n+z_up,:)
593          ENDIF
594
595          IF (ii_glo==1 .AND. jj_glo==jjm_glo) THEN
596            le(n+u_lup)=0
597            xyz_v(n+z_up,:)=xyz_v(n+z_lup,:)
598          ENDIF
599         
600        ENDDO
601      ENDDO
602
603      DO j=jj_begin-1,jj_end+1
604        DO i=ii_begin-1,ii_end+1
605          n=(j-1)*iim+i
606          xyz_i(n,:)=xyz_i(n,:) * radius
607          xyz_v(n+z_up,:)=xyz_v(n+z_up,:) * radius
608          xyz_v(n+z_down,:)=xyz_v(n+z_down,:) *radius
609          de(n+u_right)=de(n+u_right) * radius
610          de(n+u_lup)=de(n+u_lup)*radius
611          de(n+u_ldown)=de(n+u_ldown)*radius
612          xyz_e(n+u_right,:)=xyz_e(n+u_right,:)*radius
613          xyz_e(n+u_lup,:)=xyz_e(n+u_lup,:)*radius
614          xyz_e(n+u_ldown,:)=xyz_e(n+u_ldown,:)*radius
615          le(n+u_right)=le(n+u_right)*radius
616          le(n+u_lup)=le(n+u_lup)*radius
617          le(n+u_ldown)=le(n+u_ldown)*radius
618          Ai(n)=Ai(n)*radius**2
619          Av(n+z_up)=Av(n+z_up)*radius**2
620          Av(n+z_down)=Av(n+z_down)*radius**2
621        ENDDO
622      ENDDO
623                 
624    ENDDO
625
626    CALL transfert_request(geom%Ai,req_i1)
627    CALL transfert_request(geom%centroid,req_i1)
628
629!    CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)
630 
631  END SUBROUTINE set_geometry
632 
633  SUBROUTINE compute_wee(n,pos,w)
634  IMPLICIT NONE
635    INTEGER,INTENT(IN) :: n
636    INTEGER,INTENT(IN) :: pos
637    REAL(rstd),INTENT(OUT) ::w(6)
638
639    REAL(rstd) :: ne_(0:5)
640    REAL(rstd) :: Riv_(6)
641    INTEGER    :: k
642   
643   
644    DO k=0,5
645      ne_(k)=ne(n,MOD(pos-1+k+6,6)+1) 
646      Riv_(k+1)=Riv(n,MOD(pos-1+k+6,6)+1)
647    ENDDO
648         
649    w(1)=-ne_(0)*ne_(1)*(Riv_(1)-0.5)
650    w(2)=-ne_(2)*(ne_(0)*Riv_(2)-w(1)*ne_(1))
651    w(3)=-ne_(3)*(ne_(0)*Riv_(3)-w(2)*ne_(2))
652    w(4)=-ne_(4)*(ne_(0)*Riv_(4)-w(3)*ne_(3))
653    w(5)=-ne_(5)*(ne_(0)*Riv_(5)-w(4)*ne_(4))
654    w(6)=ne_(0)*ne_(5)*(Riv_(6)-0.5)
655   
656!    IF ( ABS(w(5)-w(6))>1e-20) PRINT *, "pb pour wee : w(5)!=w(6)",sum(Riv_(:))
657
658   END SUBROUTINE compute_wee
659           
660
661 
662  SUBROUTINE compute_geometry
663  IMPLICIT NONE
664    CALL allocate_geometry
665    CALL set_geometry
666   
667  END SUBROUTINE compute_geometry
668 
669END MODULE geometry
Note: See TracBrowser for help on using the repository browser.