Changeset 863 for codes/icosagcm/devel/src/sphere/geometry.f90
- Timestamp:
- 05/09/19 01:37:55 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.