Changeset 863 for codes/icosagcm/devel/src
- Timestamp:
- 05/09/19 01:37:55 (5 years ago)
- Location:
- codes/icosagcm/devel/src
- Files:
-
- 3 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
codes/icosagcm/devel/src/base/init_grid_param.f90
r830 r863 7 7 USE grid_param 8 8 USE ioipsl, ONLY : getin 9 USE init_unstructured_mod, ONLY : read_dump_partition9 USE init_unstructured_mod, ONLY : open_local_mesh_file 10 10 LOGICAL :: is_mpi_root 11 11 CHARACTER(len=255) :: grid_type_var … … 24 24 ! is_omp_level_master=.TRUE. 25 25 ! omp_level_size=1 26 CALL read_dump_partition26 CALL open_local_mesh_file 27 27 IF (is_mpi_root) PRINT *,'DYNAMICO mesh is unstructured/LAM.' 28 28 CASE DEFAULT -
codes/icosagcm/devel/src/sphere/compute_geometry.f90
r860 r863 1 MODULE geometry2 USE field_mod1 MODULE compute_geometry_mod 2 USE geometry 3 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 4 PRIVATE 5 SAVE 6 7 PUBLIC :: compute_geometry 100 8 101 9 CONTAINS 102 103 SUBROUTINE allocate_geometry104 USE field_mod105 IMPLICIT NONE106 107 CALL allocate_field(geom%Ai,field_t,type_real,name='Ai')108 109 CALL allocate_field(geom%xyz_i,field_t,type_real,3)110 CALL allocate_field(geom%lon_i,field_t,type_real)111 CALL allocate_field(geom%lat_i,field_t,type_real)112 CALL allocate_field(geom%elon_i,field_t,type_real,3)113 CALL allocate_field(geom%elat_i,field_t,type_real,3)114 CALL allocate_field(geom%centroid,field_t,type_real,3)115 116 CALL allocate_field(geom%xyz_e,field_u,type_real,3)117 CALL allocate_field(geom%lon_e,field_u,type_real)118 CALL allocate_field(geom%lat_e,field_u,type_real)119 CALL allocate_field(geom%elon_e,field_u,type_real,3)120 CALL allocate_field(geom%elat_e,field_u,type_real,3)121 CALL allocate_field(geom%ep_e,field_u,type_real,3)122 CALL allocate_field(geom%et_e,field_u,type_real,3)123 124 CALL allocate_field(geom%xyz_v,field_z,type_real,3)125 CALL allocate_field(geom%de,field_u,type_real)126 CALL allocate_field(geom%le,field_u,type_real)127 CALL allocate_field(geom%le_de,field_u,type_real)128 CALL allocate_field(geom%bi,field_t,type_real)129 CALL allocate_field(geom%Av,field_z,type_real)130 CALL allocate_field(geom%S1,field_t,type_real,6)131 CALL allocate_field(geom%S2,field_t,type_real,6)132 CALL allocate_field(geom%Riv,field_t,type_real,6)133 CALL allocate_field(geom%Riv2,field_t,type_real,6)134 CALL allocate_field(geom%ne,field_t,type_integer,6)135 CALL allocate_field(geom%Wee,field_u,type_real,5,2)136 CALL allocate_field(geom%bi,field_t,type_real)137 CALL allocate_field(geom%fv,field_z,type_real)138 139 END SUBROUTINE allocate_geometry140 141 142 SUBROUTINE swap_geometry(ind)143 USE field_mod144 USE domain_mod, ONLY : swap_needed145 IMPLICIT NONE146 INTEGER,INTENT(IN) :: ind147 148 IF(.NOT. swap_needed) RETURN149 150 !!$OMP MASTER151 Ai=geom%Ai(ind)152 xyz_i=geom%xyz_i(ind)153 centroid=geom%centroid(ind)154 xyz_e=geom%xyz_e(ind)155 ep_e=geom%ep_e(ind)156 et_e=geom%et_e(ind)157 lon_i=geom%lon_i(ind)158 lat_i=geom%lat_i(ind)159 lon_e=geom%lon_e(ind)160 lat_e=geom%lat_e(ind)161 elon_i=geom%elon_i(ind)162 elat_i=geom%elat_i(ind)163 elon_e=geom%elon_e(ind)164 elat_e=geom%elat_e(ind)165 xyz_v=geom%xyz_v(ind)166 de=geom%de(ind)167 le=geom%le(ind)168 le_de=geom%le_de(ind)169 Av=geom%Av(ind)170 S1=geom%S1(ind)171 S2=geom%S2(ind)172 Riv=geom%Riv(ind)173 Riv2=geom%Riv2(ind)174 ne=geom%ne(ind)175 Wee=geom%Wee(ind)176 bi=geom%bi(ind)177 fv=geom%fv(ind)178 !!$OMP END MASTER179 !!$OMP BARRIER180 END SUBROUTINE swap_geometry181 10 182 11 SUBROUTINE update_circumcenters … … 729 558 730 559 SUBROUTINE compute_geometry 731 IMPLICIT NONE 560 USE grid_param 561 USE domain_mod, ONLY : swap_needed 562 USE init_unstructured_mod, ONLY : read_local_mesh 563 IMPLICIT NONE 564 732 565 CALL allocate_geometry 733 CALL set_geometry 734 566 567 SELECT CASE(grid_type) 568 CASE(grid_ico) 569 CALL set_geometry 570 CASE(grid_unst) 571 swap_needed=.FALSE. 572 CALL read_local_mesh 573 CASE DEFAULT 574 STOP 'Invalid value of grid_type encountered in compute_geometry' 575 END SELECT 735 576 END SUBROUTINE compute_geometry 736 577 737 END MODULE geometry578 END MODULE compute_geometry_mod -
codes/icosagcm/devel/src/sphere/geometry.f90
r856 r863 180 180 END SUBROUTINE swap_geometry 181 181 182 SUBROUTINE update_circumcenters183 USE domain_mod184 USE dimensions185 USE spherical_geom_mod186 USE vector187 USE transfert_mod188 USE omp_para189 190 IMPLICIT NONE191 REAL(rstd) :: x1(3),x2(3)192 REAL(rstd) :: vect(3,6)193 REAL(rstd) :: centr(3)194 INTEGER :: ind,i,j,n,k195 TYPE(t_message),SAVE :: message0, message1196 LOGICAL, SAVE :: first=.TRUE.197 !$OMP THREADPRIVATE(first)198 199 IF (first) THEN200 CALL init_message(geom%xyz_i, req_i0 ,message0)201 CALL init_message(geom%xyz_i, req_i1 ,message1)202 first=.FALSE.203 ENDIF204 205 CALL transfert_message(geom%xyz_i,message0)206 CALL transfert_message(geom%xyz_i,message1)207 208 DO ind=1,ndomain209 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE210 CALL swap_dimensions(ind)211 CALL swap_geometry(ind)212 DO j=jj_begin,jj_end213 DO i=ii_begin,ii_end214 n=(j-1)*iim+i215 DO k=0,5216 x1(:) = xyz_i(n+t_pos(k+1),:)217 x2(:) = xyz_i(n+t_pos(MOD(k+1,6)+1),:)218 if (norm(x1-x2)<1e-16) x2(:) = xyz_i(n+t_pos(MOD(k+2,6)+1),:)219 CALL circumcenter(xyz_i(n,:), x1, x2, xyz_v(n+z_pos(k+1),:))220 ENDDO221 ENDDO222 ENDDO223 ENDDO224 225 END SUBROUTINE update_circumcenters226 227 SUBROUTINE remap_schmidt_loc228 USE spherical_geom_mod229 USE getin_mod230 USE omp_para231 USE domain_mod232 USE dimensions233 IMPLICIT NONE234 INTEGER :: ind,i,j,n235 REAL(rstd) :: schmidt_factor, schmidt_lon, schmidt_lat236 237 ! Schmidt transform parameters238 schmidt_factor = 1.239 CALL getin('schmidt_factor', schmidt_factor)240 schmidt_factor = schmidt_factor**2.241 242 schmidt_lon = 0.243 CALL getin('schmidt_lon', schmidt_lon)244 schmidt_lon = schmidt_lon * pi/180.245 246 schmidt_lat = 45.247 CALL getin('schmidt_lat', schmidt_lat)248 schmidt_lat = schmidt_lat * pi/180.249 250 DO ind=1,ndomain251 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE252 CALL swap_dimensions(ind)253 CALL swap_geometry(ind)254 DO j=jj_begin,jj_end255 DO i=ii_begin,ii_end256 n=(j-1)*iim+i257 CALL schmidt_transform(xyz_i(n,:), schmidt_factor, schmidt_lon, schmidt_lat)258 ENDDO259 ENDDO260 ENDDO261 END SUBROUTINE remap_schmidt_loc262 263 SUBROUTINE optimize_geometry264 USE metric265 USE spherical_geom_mod266 USE domain_mod267 USE dimensions268 USE transfert_mod269 USE vector270 USE getin_mod271 USE omp_para272 IMPLICIT NONE273 INTEGER :: nb_it=0274 TYPE(t_domain),POINTER :: d275 INTEGER :: ind,it,i,j,n,k276 REAL(rstd) :: x1(3),x2(3)277 REAL(rstd) :: vect(3,6)278 REAL(rstd) :: centr(3)279 REAL(rstd) :: sum280 LOGICAL :: check281 282 283 CALL getin('optim_it',nb_it)284 285 DO ind=1,ndomain286 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE287 d=>domain(ind)288 CALL swap_dimensions(ind)289 CALL swap_geometry(ind)290 xyz_i(:,1) = 0 ; xyz_i(:,2) = 0 ; xyz_i(:,3) = 1291 292 DO j=jj_begin,jj_end293 DO i=ii_begin,ii_end294 n=(j-1)*iim+i295 xyz_i(n,:)=d%xyz(:,i,j)296 ENDDO297 ENDDO298 ENDDO299 300 CALL update_circumcenters301 302 DO ind=1,ndomain303 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE304 d=>domain(ind)305 CALL swap_dimensions(ind)306 CALL swap_geometry(ind)307 DO j=jj_begin,jj_end308 DO i=ii_begin,ii_end309 n=(j-1)*iim+i310 DO k=0,5311 x1(:) = xyz_v(n+z_pos(k+1),:)312 x2(:) = d%vertex(:,k,i,j)313 IF (norm(x1-x2)>1e-10) THEN314 PRINT*,"vertex diff ",ind,i,j,k315 PRINT*,x1316 PRINT*,x2317 ENDIF318 ENDDO319 ENDDO320 ENDDO321 ENDDO322 323 324 DO it=1,nb_it325 IF (MOD(it,100)==0) THEN326 check=is_master327 ELSE328 check=.FALSE.329 ENDIF330 331 sum=0332 DO ind=1,ndomain333 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE334 CALL swap_dimensions(ind)335 CALL swap_geometry(ind)336 DO j=jj_begin,jj_end337 DO i=ii_begin,ii_end338 n=(j-1)*iim+i339 vect(:,1)=xyz_v(n+z_rup,:)340 vect(:,2)=xyz_v(n+z_up,:)341 vect(:,3)=xyz_v(n+z_lup,:)342 vect(:,4)=xyz_v(n+z_ldown,:)343 vect(:,5)=xyz_v(n+z_down,:)344 vect(:,6)=xyz_v(n+z_rdown,:)345 CALL compute_centroid(vect,6,centr)346 IF (check) THEN347 sum=MAX(sum,norm(xyz_i(n,:)-centr(:)))348 ENDIF349 xyz_i(n,:)=centr(:)350 ENDDO351 ENDDO352 353 ENDDO354 355 IF (check) PRINT *,"it = ",it," diff centroid circumcenter ",sum356 357 CALL update_circumcenters358 359 ENDDO360 361 END SUBROUTINE optimize_geometry362 363 SUBROUTINE update_domain364 ! copy position of generators and vertices back into domain(:)%xyz/vertex365 ! so that XIOS/create_header_gen gets the right values366 USE omp_para367 USE dimensions368 USE domain_mod369 USE transfert_mpi_mod370 371 INTEGER :: ind,i,j,k,n372 TYPE(t_domain),POINTER :: d373 TYPE(t_field),POINTER,SAVE :: xyz_glo(:), xyz_loc(:), vertex_glo(:), vertex_loc(:)374 REAL(rstd), POINTER :: xyz(:,:), vertex(:,:)375 376 CALL allocate_field(xyz_loc, field_t, type_real, 3)377 CALL allocate_field(vertex_loc, field_z, type_real, 3)378 379 DO ind=1,ndomain380 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE381 CALL swap_dimensions(ind)382 CALL swap_geometry(ind)383 xyz = xyz_loc(ind)384 xyz(:,:) = xyz_i(:,:)385 vertex = vertex_loc(ind)386 vertex(:,:) = xyz_v(:,:)387 END DO388 389 !$OMP BARRIER390 !$OMP MASTER391 CALL allocate_field_glo(xyz_glo, field_t, type_real, 3)392 CALL allocate_field_glo(vertex_glo, field_z, type_real, 3)393 394 CALL gather_field(xyz_loc, xyz_glo)395 CALL gather_field(vertex_loc, vertex_glo)396 397 CALL bcast_field(xyz_glo)398 CALL bcast_field(vertex_glo)399 400 DO ind=1,ndomain_glo401 d=>domain_glo(ind)402 xyz = xyz_glo(ind)403 vertex = vertex_glo(ind)404 DO j=d%jj_begin,d%jj_end405 DO i=d%ii_begin,d%ii_end406 n=(j-1)*d%iim+i407 d%xyz(:,i,j)=xyz(n,:)408 DO k=0,5409 d%vertex(:,k,i,j) = vertex(n+d%z_pos(k+1),:)410 END DO411 END DO412 END DO413 END DO414 415 CALL deallocate_field_glo(vertex_glo)416 CALL deallocate_field_glo(xyz_glo)417 !$OMP END MASTER418 !$OMP BARRIER419 CALL deallocate_field(vertex_loc)420 CALL deallocate_field(xyz_loc)421 422 END SUBROUTINE update_domain423 424 SUBROUTINE set_geometry425 USE metric426 USE vector427 USE spherical_geom_mod428 USE domain_mod429 USE dimensions430 USE transfert_mod431 USE getin_mod432 USE omp_para433 IMPLICIT NONE434 435 REAL(rstd) :: surf(6)436 REAL(rstd) :: surf_v(6)437 REAL(rstd) :: vect(3,6)438 REAL(rstd) :: centr(3)439 REAL(rstd) :: vet(3),vep(3), vertex(3)440 INTEGER :: ind,i,j,k,n441 TYPE(t_domain),POINTER :: d442 REAL(rstd) :: S12443 REAL(rstd) :: w(6)444 REAL(rstd) :: lon,lat445 INTEGER :: ii_glo,jj_glo446 REAL(rstd) :: S447 448 449 CALL optimize_geometry450 CALL remap_schmidt_loc451 CALL update_circumcenters452 ! copy position of generators and vertices back into domain(:)%xyz/vertex453 ! so that XIOS gets the right values454 CALL update_domain455 456 DO ind=1,ndomain457 IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master ) CYCLE458 d=>domain(ind)459 CALL swap_dimensions(ind)460 CALL swap_geometry(ind)461 lon_i(:)=0 ; lat_i(:)=0462 lon_e(:)=0 ; lat_e(:)=0463 DO j=jj_begin-1,jj_end+1464 DO i=ii_begin-1,ii_end+1465 n=(j-1)*iim+i466 467 DO k=0,5468 ne(n,k+1)=d%ne(k,i,j)469 ENDDO470 471 vect(:,1)=xyz_v(n+z_rup,:)472 vect(:,2)=xyz_v(n+z_up,:)473 vect(:,3)=xyz_v(n+z_lup,:)474 vect(:,4)=xyz_v(n+z_ldown,:)475 vect(:,5)=xyz_v(n+z_down,:)476 vect(:,6)=xyz_v(n+z_rdown,:)477 CALL compute_centroid(vect,6,centr)478 centroid(n,:)=centr(:)479 480 481 CALL xyz2lonlat(xyz_v(n+z_up,:),lon,lat)482 fv(n+z_up)=2*sin(lat)*omega483 CALL xyz2lonlat(xyz_v(n+z_down,:),lon,lat)484 fv(n+z_down)=2*sin(lat)*omega485 486 bi(n)=0.487 488 CALL dist_cart(xyz_i(n,:),xyz_i(n+t_right,:),de(n+u_right))489 CALL dist_cart(xyz_i(n,:),xyz_i(n+t_lup,:),de(n+u_lup))490 CALL dist_cart(xyz_i(n,:),xyz_i(n+t_ldown,:),de(n+u_ldown))491 492 CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_right,:),0.5,xyz_e(n+u_right,:))493 CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_lup,:),0.5,xyz_e(n+u_lup,:))494 CALL div_arc_bis(xyz_i(n,:),xyz_i(n+t_ldown,:),0.5,xyz_e(n+u_ldown,:))495 496 CALL dist_cart(xyz_v(n+z_rdown,:), xyz_v(n+z_rup,:),le(n+u_right))497 CALL dist_cart(xyz_v(n+z_up,:), xyz_v(n+z_lup,:),le(n+u_lup))498 CALL dist_cart(xyz_v(n+z_ldown,:), xyz_v(n+z_down,:),le(n+u_ldown))499 500 le_de(n+u_right)=le(n+u_right)/de(n+u_right) ! NaN possible but should be harmless501 le_de(n+u_lup) =le(n+u_lup) /de(n+u_lup)502 le_de(n+u_ldown)=le(n+u_ldown)/de(n+u_ldown)503 504 Ai(n)=0505 DO k=0,5506 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))507 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))508 Ai(n)=Ai(n)+surf(k+1)509 IF (i==ii_end .AND. j==jj_begin) THEN510 IF (Ai(n)<1e20) THEN511 ELSE512 PRINT *,"PB !!",Ai(n),k,surf(k+1)513 PRINT*,xyz_i(n,:),xyz_v(n+z_pos(MOD((k-1+6),6)+1),:),xyz_v(n+z_pos(k+1),:)514 ENDIF515 ENDIF516 ENDDO517 518 ! Sign convention : Ringler et al., JCP 2010, eq. 21 p. 3071519 ! Normal component is along outgoing normal vector if ne=1520 521 CALL cross_product2(xyz_v(n+z_rdown,:),xyz_v(n+z_rup,:),vep)522 IF (norm(vep)>1e-30) THEN523 vep(:)=vep(:)/norm(vep) ! Inward normal vector524 CALL cross_product2(vep,xyz_e(n+u_right,:),vet) ! Counter-clockwise tangent vector525 vet(:)=vet(:)/norm(vet)526 ep_e(n+u_right,:)=-vep(:)*ne(n,right)527 et_e(n+u_right,:)=vet(:)*ne(n,right)528 ENDIF529 530 CALL cross_product2(xyz_v(n+z_up,:),xyz_v(n+z_lup,:),vep)531 IF (norm(vep)>1e-30) THEN532 vep(:)=vep(:)/norm(vep)533 CALL cross_product2(vep,xyz_e(n+u_lup,:),vet)534 vet(:)=vet(:)/norm(vet)535 ep_e(n+u_lup,:)=-vep(:)*ne(n,lup)536 et_e(n+u_lup,:)=vet(:)*ne(n,lup)537 ENDIF538 539 CALL cross_product2(xyz_v(n+z_ldown,:),xyz_v(n+z_down,:),vep)540 IF (norm(vep)>1e-30) THEN541 vep(:)=vep(:)/norm(vep)542 CALL cross_product2(vep,xyz_e(n+u_ldown,:),vet)543 vet(:)=vet(:)/norm(vet)544 ep_e(n+u_ldown,:)=-vep(:)*ne(n,ldown)545 et_e(n+u_ldown,:)=vet(:)*ne(n,ldown)546 ENDIF547 548 CALL xyz2lonlat(xyz_i(n,:),lon,lat)549 lon_i(n)=lon550 lat_i(n)=lat551 elon_i(n,1) = -sin(lon)552 elon_i(n,2) = cos(lon)553 elon_i(n,3) = 0554 elat_i(n,1) = -cos(lon)*sin(lat)555 elat_i(n,2) = -sin(lon)*sin(lat)556 elat_i(n,3) = cos(lat)557 558 559 CALL xyz2lonlat(xyz_e(n+u_right,:),lon,lat)560 lon_e(n+u_right)=lon561 lat_e(n+u_right)=lat562 elon_e(n+u_right,1) = -sin(lon)563 elon_e(n+u_right,2) = cos(lon)564 elon_e(n+u_right,3) = 0565 elat_e(n+u_right,1) = -cos(lon)*sin(lat)566 elat_e(n+u_right,2) = -sin(lon)*sin(lat)567 elat_e(n+u_right,3) = cos(lat)568 569 CALL xyz2lonlat(xyz_e(n+u_lup,:),lon,lat)570 lon_e(n+u_lup)=lon571 lat_e(n+u_lup)=lat572 elon_e(n+u_lup,1) = -sin(lon)573 elon_e(n+u_lup,2) = cos(lon)574 elon_e(n+u_lup,3) = 0575 elat_e(n+u_lup,1) = -cos(lon)*sin(lat)576 elat_e(n+u_lup,2) = -sin(lon)*sin(lat)577 elat_e(n+u_lup,3) = cos(lat)578 579 CALL xyz2lonlat(xyz_e(n+u_ldown,:),lon,lat)580 lon_e(n+u_ldown)=lon581 lat_e(n+u_ldown)=lat582 elon_e(n+u_ldown,1) = -sin(lon)583 elon_e(n+u_ldown,2) = cos(lon)584 elon_e(n+u_ldown,3) = 0585 elat_e(n+u_ldown,1) = -cos(lon)*sin(lat)586 elat_e(n+u_ldown,2) = -sin(lon)*sin(lat)587 elat_e(n+u_ldown,3) = cos(lat)588 589 590 DO k=0,5591 CALL surf_triangle(xyz_i(n,:), xyz_v(n+z_pos(k+1),:), xyz_i(n+t_pos(k+1),:),S1(n,k+1) )592 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(n,k+1) )593 S12 = .5*(S1(n,k+1)+S2(n,k+1))594 Riv(n,k+1)=S12/Ai(n)595 Riv2(n,k+1)=S12/surf_v(k+1)596 ENDDO597 598 DO k=1,6599 IF (ABS(surf_v(k))<1e-30) THEN600 Riv(n,k)=0.601 ENDIF602 ENDDO603 604 Av(n+z_up)=surf_v(vup)+1e-100605 Av(n+z_down)=surf_v(vdown)+1e-100606 607 ENDDO608 ENDDO609 610 DO j=jj_begin,jj_end611 DO i=ii_begin,ii_end612 n=(j-1)*iim+i613 614 CALL compute_wee(n,right,w)615 Wee(n+u_right,:,1)=w(1:5)616 617 CALL compute_wee(n+t_right,left,w)618 Wee(n+u_right,:,2)=w(1:5)619 620 621 CALL compute_wee(n,lup,w)622 Wee(n+u_lup,:,1)=w(1:5)623 624 CALL compute_wee(n+t_lup,rdown,w)625 Wee(n+u_lup,:,2)=w(1:5)626 627 628 CALL compute_wee(n,ldown,w)629 Wee(n+u_ldown,:,1)=w(1:5)630 631 CALL compute_wee(n+t_ldown,rup,w)632 Wee(n+u_ldown,:,2)=w(1:5)633 634 ENDDO635 ENDDO636 637 DO j=jj_begin,jj_end638 DO i=ii_begin,ii_end639 n=(j-1)*iim+i640 ii_glo=d%ii_begin_glo-d%ii_begin+i641 jj_glo=d%jj_begin_glo-d%jj_begin+j642 643 IF (ii_glo==1 .AND. jj_glo==1) THEN644 le(n+u_ldown)=0645 le_de(n+u_ldown)=0646 xyz_v(n+z_ldown,:)=xyz_v(n+z_down,:)647 648 ENDIF649 650 IF (ii_glo==iim_glo .AND. jj_glo==1) THEN651 le(n+u_right)=0652 le_de(n+u_right)=0653 xyz_v(n+z_rdown,:)=xyz_v(n+z_rup,:)654 ENDIF655 656 IF (ii_glo==iim_glo .AND. jj_glo==jjm_glo) THEN657 le(n+u_rup)=0658 le_de(n+u_rup)=0659 xyz_v(n+z_rup,:)=xyz_v(n+z_up,:)660 ENDIF661 662 IF (ii_glo==1 .AND. jj_glo==jjm_glo) THEN663 le(n+u_lup)=0664 le_de(n+u_lup)=0665 xyz_v(n+z_up,:)=xyz_v(n+z_lup,:)666 ENDIF667 668 ENDDO669 ENDDO670 671 DO j=jj_begin-1,jj_end+1672 DO i=ii_begin-1,ii_end+1673 n=(j-1)*iim+i674 xyz_i(n,:)=xyz_i(n,:) * radius675 xyz_v(n+z_up,:)=xyz_v(n+z_up,:) * radius676 xyz_v(n+z_down,:)=xyz_v(n+z_down,:) *radius677 de(n+u_right)=de(n+u_right) * radius678 de(n+u_lup)=de(n+u_lup)*radius679 de(n+u_ldown)=de(n+u_ldown)*radius680 xyz_e(n+u_right,:)=xyz_e(n+u_right,:)*radius681 xyz_e(n+u_lup,:)=xyz_e(n+u_lup,:)*radius682 xyz_e(n+u_ldown,:)=xyz_e(n+u_ldown,:)*radius683 le(n+u_right)=le(n+u_right)*radius684 le(n+u_lup)=le(n+u_lup)*radius685 le(n+u_ldown)=le(n+u_ldown)*radius686 Ai(n)=Ai(n)*radius**2687 Av(n+z_up)=Av(n+z_up)*radius**2688 Av(n+z_down)=Av(n+z_down)*radius**2689 ENDDO690 ENDDO691 692 ENDDO693 694 CALL transfert_request(geom%Ai,req_i1)695 CALL transfert_request(geom%centroid,req_i1)696 697 ! CALL surf_triangle(d%xyz(:,ii_begin,jj_begin),d%xyz(:,ii_begin,jj_end),d%xyz(:,ii_end,jj_begin),S)698 699 END SUBROUTINE set_geometry700 701 SUBROUTINE compute_wee(n,pos,w)702 IMPLICIT NONE703 INTEGER,INTENT(IN) :: n704 INTEGER,INTENT(IN) :: pos705 REAL(rstd),INTENT(OUT) ::w(6)706 707 REAL(rstd) :: ne_(0:5)708 REAL(rstd) :: Riv_(6)709 INTEGER :: k710 711 712 DO k=0,5713 ne_(k)=ne(n,MOD(pos-1+k+6,6)+1)714 Riv_(k+1)=Riv(n,MOD(pos-1+k+6,6)+1)715 ENDDO716 717 w(1)=-ne_(0)*ne_(1)*(Riv_(1)-0.5)718 w(2)=-ne_(2)*(ne_(0)*Riv_(2)-w(1)*ne_(1))719 w(3)=-ne_(3)*(ne_(0)*Riv_(3)-w(2)*ne_(2))720 w(4)=-ne_(4)*(ne_(0)*Riv_(4)-w(3)*ne_(3))721 w(5)=-ne_(5)*(ne_(0)*Riv_(5)-w(4)*ne_(4))722 w(6)=ne_(0)*ne_(5)*(Riv_(6)-0.5)723 724 ! IF ( ABS(w(5)-w(6))>1e-20) PRINT *, "pb pour wee : w(5)!=w(6)",sum(Riv_(:))725 726 END SUBROUTINE compute_wee727 728 729 730 SUBROUTINE compute_geometry731 IMPLICIT NONE732 CALL allocate_geometry733 CALL set_geometry734 735 END SUBROUTINE compute_geometry736 737 182 END MODULE geometry -
codes/icosagcm/devel/src/unstructured/init_unstructured.f90
r856 r863 2 2 USE mpipara, ONLY : is_mpi_master 3 3 USE data_unstructured_mod 4 USE geometry, ONLY : de 4 5 IMPLICIT NONE 5 6 SAVE … … 10 11 INTEGER, ALLOCATABLE :: Idata_read1(:),Idata_read2(:,:),Idata_read3(:,:,:) 11 12 12 PUBLIC :: read_dump_partition 13 CHARACTER(LEN=*),PARAMETER :: meshfile="input/mesh_information.nc" 14 INTEGER :: id_nc ! NetCDF id of mesh file open by open_local_mesh_file 15 16 PUBLIC :: open_local_mesh_file, read_local_mesh 13 17 14 18 CONTAINS … … 92 96 93 97 94 SUBROUTINE read_dump_partition 95 use netcdf_mod 96 USE ioipsl 97 USE field_mod 98 USE geometry, ONLY : lon_i, lat_i, lon_e, lat_e, ep_e 99 IMPLICIT NONE 100 101 !!-------------Declare local variables------------------- 102 CHARACTER(LEN=*),PARAMETER :: f="input/mesh_information.nc" 103 INTEGER :: id_nc, ierr, status, descriptionLength, ij 98 SUBROUTINE open_local_mesh_file 99 USE netcdf_mod 104 100 CHARACTER(LEN= 80) :: description 105 REAL(rstd), ALLOCATABLE :: angle_e(:) 106 REAL(rstd) :: clon, slon, clat, slat, & ! COS/SIN of lon/lat 107 elon(3), elat(3) ! lon/lat unit vectors 108 print *,"------------------ READING FILE " , f, "----------------------- " 101 INTEGER :: ierr, status, descriptionLength 102 103 PRINT *,"------------------ READING FILE " , meshfile, "----------------------- " 109 104 !open and read the input file 110 ierr = nf90_open( f, NF90_NOWRITE, id_nc)105 ierr = nf90_open(meshfile, NF90_NOWRITE, id_nc) 111 106 if (ierr /= nf90_noerr) then 112 107 print *, trim(nf90_strerror(ierr)) … … 114 109 end if 115 110 116 status = nf90_inquire_attribute(id_nc, nf90_global, "description", len 111 status = nf90_inquire_attribute(id_nc, nf90_global, "description", len=descriptionLength) 117 112 IF(status /= 0 .or. len(description) < descriptionLength) THEN 118 113 print *, "Not enough space to put NetCDF attribute values." … … 124 119 print *,"Data file description :",description 125 120 121 CALL read_from_gridfile(id_nc, 'integer', 'primal_deg') 122 primal_num = SIZE(Idata_read1) 123 CALL read_from_gridfile(id_nc, 'integer', 'dual_deg') 124 dual_num = SIZE(Idata_read1) 125 CALL read_from_gridfile(id_nc, 'integer', 'trisk_deg') 126 edge_num = SIZE(Idata_read1) 127 END SUBROUTINE open_local_mesh_file 128 129 130 SUBROUTINE read_local_mesh 131 USE field_mod 132 USE geometry, ONLY : lon_i, lat_i, lon_e, lat_e, ep_e 133 IMPLICIT NONE 134 INTEGER :: ij 135 REAL(rstd), ALLOCATABLE :: angle_e(:) 136 REAL(rstd) :: clon, slon, clat, slat, & ! COS/SIN of lon/lat 137 elon(3), elat(3) ! lon/lat unit vectors 138 126 139 !status = nf90_get_att(id_nc, nf90_global, "primal_num", primal_num) 127 140 !status = nf90_get_att(id_nc, nf90_global, "dual_num", dual_num) … … 161 174 CALL read_from_gridfile(id_nc, 'float', 'le_de') 162 175 ALLOCATE(le_de, source = Ddata_read1) 176 CALL read_from_gridfile(id_nc, 'float', 'le') 177 ALLOCATE(le, source = Ddata_read1) 178 CALL read_from_gridfile(id_nc, 'float', 'de') 179 ALLOCATE(de, source = Ddata_read1) 163 180 CALL read_from_gridfile(id_nc, 'float', 'Riv2') 164 181 ALLOCATE(Riv2, source = Ddata_read2) … … 204 221 205 222 DEALLOCATE(angle_e) 206 END SUBROUTINE read_dump_partition 207 208 SUBROUTINE init_grid_type 209 USE grid_param, ONLY : grid_type, grid_unst, grid_ico 210 USE getin_mod, ONLY : getin 211 CHARACTER(len=255) :: grid_type_var 212 grid_type_var='icosahedral' 213 CALL getin("grid_type",grid_type_var) 214 SELECT CASE(grid_type_var) 215 CASE('unstructured') 216 grid_type = grid_unst 217 ! is_omp_level_master=.TRUE. 218 ! omp_level_size=1 219 CALL read_dump_partition 220 IF (is_mpi_master) PRINT *,'Using unstructured grid type' 221 CASE DEFAULT 222 grid_type = grid_ico 223 IF (is_mpi_master) PRINT *,'Using default grid type' 224 END SELECT 225 END SUBROUTINE init_grid_type 223 END SUBROUTINE read_local_mesh 226 224 227 225 END MODULE init_unstructured_mod
Note: See TracChangeset
for help on using the changeset viewer.