Changeset 8365
- Timestamp:
- 2024-01-08T15:20:57+01:00 (4 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flow.f90
r8300 r8365 24 24 25 25 !! PARAMETERS 26 REAL(r_std), SAVE :: fast_tcst = 3.0 !! Property of the fast reservoir (day/m)26 REAL(r_std), SAVE :: fast_tcst = 9e-3 !! Property of the fast reservoir (day/km) 27 27 !$OMP THREADPRIVATE(fast_tcst) 28 REAL(r_std), SAVE :: slow_tcst = 25.0 !! Property of the slow reservoir (day/m)28 REAL(r_std), SAVE :: slow_tcst = 1.2e-2 !! Property of the slow reservoir (day/km) 29 29 !$OMP THREADPRIVATE(slow_tcst) 30 REAL(r_std), SAVE :: stream_tcst = 0.24 !! Property of the stream reservoir (day/m)30 REAL(r_std), SAVE :: stream_tcst = 3e-5 !! Property of the stream reservoir (day/km) 31 31 !$OMP THREADPRIVATE(stream_tcst) 32 32 … … 41 41 !$OMP THREADPRIVATE(fast_reservoir_r) 42 42 43 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET:: slow_reservoir_r(:) !! Water amount in the slow reservoir (kg) - (local routing grid)43 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET:: slow_reservoir_r(:) !! Water amount in the slow reservoir (kg) - (local routing grid) 44 44 !$OMP THREADPRIVATE(slow_reservoir_r) 45 45 … … 47 47 !$OMP THREADPRIVATE(stream_reservoir_r) 48 48 49 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: riverflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid)49 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: riverflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid) 50 50 !$OMP THREADPRIVATE(riverflow_mean) 51 51 … … 53 53 !$OMP THREADPRIVATE(coastalflow_mean) 54 54 55 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: lakeinflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid)55 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: lakeinflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid) 56 56 !$OMP THREADPRIVATE(lakeinflow_mean) 57 57 … … 61 61 REAL(r_std),SAVE,ALLOCATABLE :: topoind_r(:) !! Topographic index of the retention time (m) index - (local routing grid) 62 62 !$OMP THREADPRIVATE(topoind_r) 63 REAL(r_std),SAVE :: topoind_factor !! conversion factor for topographic index (merit : m, standard : km), topo index must be in m 64 !$OMP THREADPRIVATE(topoind_factor) 63 65 INTEGER,SAVE,ALLOCATABLE :: route_flow_rp1(:) !! flow index from cell to neighboring cell following the trip direction - (local routing grid + halo) 64 66 !$OMP THREADPRIVATE(route_flow_rp1) … … 69 71 LOGICAL,SAVE,ALLOCATABLE :: is_riverflow_r(:) !! is river flow point - (local routing grid) 70 72 !$OMP THREADPRIVATE(is_riverflow_r) 71 LOGICAL,SAVE,ALLOCATABLE :: is_coastline(:) !! is coastline point - (local native grid)73 LOGICAL,SAVE,ALLOCATABLE :: is_coastline(:) !! is coastline point - (local native grid) 72 74 !$OMP THREADPRIVATE(is_coastline) 73 75 LOGICAL,SAVE,ALLOCATABLE :: is_streamflow_r(:) !! is stream flow point - (local routing grid) … … 85 87 REAL(r_std),SAVE,ALLOCATABLE :: routing_weight(:) !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid) 86 88 !$OMP THREADPRIVATE(routing_weight) 87 REAL(r_std),SAVE,ALLOCATABLE :: routing_weight_in(:) 89 REAL(r_std),SAVE,ALLOCATABLE :: routing_weight_in(:) !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid) 88 90 !$OMP THREADPRIVATE(routing_weight_in) 89 91 REAL(r_std),SAVE,ALLOCATABLE :: unrouted_weight(:) !! Weight to transform runoff and drainage flux to water quantity in a conservative way (local native grid -> local routing grid) … … 253 255 INTEGER ,INTENT(IN) :: kjit 254 256 INTEGER ,INTENT(IN) :: rest_id 255 INTEGER, INTENT(IN) :: nbpt_ !! nb points orchidee grid257 INTEGER, INTENT(IN) :: nbpt_ !! nb points orchidee grid 256 258 REAL(r_std), INTENT(IN) :: dt_routing 257 REAL(r_std), INTENT(IN) :: contfrac(nbpt_) !! fraction of land258 INTEGER, INTENT(OUT) :: nbpt_r_ !! nb points routing native grid259 REAL(r_std), INTENT(IN) :: contfrac(nbpt_) !! fraction of land 260 INTEGER, INTENT(OUT) :: nbpt_r_ !! nb points routing native grid 259 261 REAL(r_std), INTENT(OUT) :: riverflow(nbpt_) !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt) 260 262 REAL(r_std), INTENT(OUT) :: coastalflow(nbpt_) !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt) … … 311 313 LOGICAL :: coastline_glo(nbp_glo) 312 314 LOGICAL :: coastline2D_glo(iim_g*jjm_g) 315 LOGICAL :: mask1d(nbpt) 316 LOGICAL :: mask1d_glo(nbp_glo) 317 LOGICAL :: mask2d_glo(iim_g*jjm_g) 313 318 INTEGER(i_std) :: i,j,ij,next_i,next_j,next_ij,m,k 314 319 REAL(r_std), PARAMETER :: epsilon=1e-5 … … 316 321 317 322 is_periodic = global 323 mask1d=.TRUE. ; 324 CALL gather(mask1d,mask1d_glo) 318 325 CALL gather(contfrac,contfrac_glo) 319 contfrac2D_glo(:)=0 320 DO i=1,nbp_glo 321 contfrac2D_glo(index_g(i)) = contfrac_glo(i) 322 ENDDO 326 323 327 324 328 IF (is_mpi_root .AND. is_omp_root) THEN 325 329 contfrac2D_glo(:)=0 330 mask2d_glo=.FALSE. 331 DO i=1,nbp_glo 332 contfrac2D_glo(index_g(i)) = contfrac_glo(i) 333 mask2d_glo(index_g(i))=mask1d_glo(i) 334 ENDDO 335 326 336 coastline2D_glo(:)=.FALSE. 327 337 DO j=1,jjm_g 328 338 DO i=1,iim_g 329 339 ij=(j-1)*iim_g+i 330 IF (contfrac2D_glo(ij)>epsilon .AND. contfrac2D_glo(ij)<1-epsilon) THEN 331 coastline2D_glo(ij)=.TRUE. 332 ELSE IF (contfrac2D_glo(ij)>=1-epsilon .AND. grid_type/=unstructured ) THEN 340 IF (grid_type==unstructured) THEN 341 IF (contfrac2D_glo(ij)>epsilon .AND. contfrac2D_glo(ij)<1-epsilon) THEN 342 coastline2D_glo(ij)=.TRUE. 343 ENDIF 344 ELSE 333 345 DO k=-1,1 334 346 DO m=-1,1 … … 358 370 CYCLE 359 371 ENDIF 360 next_ij = (next_j-1)*iim_g+next_i 361 IF (contfrac2D_glo(next_ij)<=epsilon) coastline2D_glo(ij)=.TRUE. 372 next_ij = (next_j-1)*iim_g+next_i 373 IF (.NOT. mask2d_glo(next_ij)) coastline2D_glo(ij)=.TRUE. 374 ! IF (contfrac2D_glo(next_ij)<=epsilon) coastline2D_glo(ij)=.TRUE. 362 375 ENDDO 363 376 ENDDO … … 376 389 END SUBROUTINE compute_coastline 377 390 378 379 380 SUBROUTINE routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r) 391 392 SUBROUTINE new_routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r) 381 393 USE xios 382 USE grid, ONLY : global394 ! USE grid, ONLY : global 383 395 IMPLICIT NONE 384 396 INCLUDE "mpif.h" … … 391 403 REAL(r_std), INTENT(INOUT) :: topoind_r(ni,nj) ! INPUT/OUTPUT : topographic index which will be modified by the routine if extended to ocean 392 404 405 REAL(r_std) :: mask(nbp_mpi) 393 406 REAL(r_std) :: contfrac_r(ni,nj) 394 407 REAL(r_std) :: mask_coast(nbp_mpi) … … 398 411 REAL(r_std) :: state_r(ni,nj) 399 412 REAL(r_std) :: state_rp1(0:ni+1,0:nj+1) 400 413 REAL(r_std) :: done_rp1(0:ni+1,0:nj+1) 414 REAL(r_std) :: trip_out_r(0:ni+1,0:nj+1) 415 401 416 INTEGER, PARAMETER :: is_ter=0 402 417 INTEGER, PARAMETER :: is_coast=1 403 418 INTEGER, PARAMETER :: is_oce=2 419 INTEGER, PARAMETER :: riverflow=99 420 INTEGER, PARAMETER :: coastalflow=98 421 INTEGER, PARAMETER :: lakeinflow=97 422 INTEGER, PARAMETER :: norouting=100 423 404 424 INTEGER :: updated 405 425 INTEGER :: it, i,j,k,m, nexti,nextj,previ,prevj,ierr 406 426 INTEGER :: ibegin, jbegin, ni_glo, nj_glo 407 427 REAL(r_std) :: lonvalue_1d(ni), latvalue_1d(nj) 408 REAL(r_std),PARAMETER :: epsilon=1e- 3428 REAL(r_std),PARAMETER :: epsilon=1e-5 409 429 410 430 TYPE(xios_duration) :: ts … … 419 439 CHARACTER(LEN=20) :: calendar_type 420 440 LOGICAL :: true=.TRUE. 421 441 LOGICAL :: global=.FALSE. 442 LOGICAL :: ok 443 444 CALL xios_get_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni_glo=ni_glo,nj_glo=nj_glo) ! get routing domain dimension 445 CALL xios_get_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! get routing domain dimension 446 447 ! manage non periodic boundary in case of LAM 448 IF (.NOT. global) THEN 449 IF (ibegin==0) THEN 450 i=1 ; 451 DO j=1,nj 452 IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow 453 IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow 454 ENDDO 455 ENDIF 456 457 IF (ibegin+ni==ni_glo) THEN 458 i=ni ; 459 DO j=1,nj 460 IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow 461 IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow 462 ENDDO 463 ENDIF 464 465 IF (jbegin==0) THEN 466 j=1 ; 467 DO i=1,ni 468 IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow 469 IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow 470 ENDDO 471 ENDIF 472 473 IF (jbegin+nj==nj_glo) THEN 474 j=nj ; 475 DO i=1,ni 476 IF (trip_r(i,j)>=1 .AND. trip_r(i,j)<=8) trip_r(i,j) = coastalflow 477 IF (trip_extended_r(i,j)>=1 .AND. trip_extended_r(i,j)<=8) trip_extended_r(i,j) = coastalflow 478 ENDDO 479 ENDIF 480 ENDIF 481 482 mask = 1 483 CALL xios_send_field("routing_contfrac", mask) 484 CALL xios_recv_field("routing_contfrac_r", contfrac_r) 485 486 CALL xios_send_field("trip_ext_r",trip_extended_r) 487 CALL xios_recv_field("trip_ext_rp1",trip_extended_rp1) 488 489 mask_coast=0 490 WHERE(coastline) mask_coast=1 491 CALL xios_send_field("mask_coastline",mask_coast) 492 CALL xios_recv_field("frac_coastline_r",frac_coast_r) 493 494 DO j=1,nj 495 DO i=1,ni 496 IF (frac_coast_r(i,j)>epsilon) THEN 497 state_r(i,j)=is_coast 498 ELSE 499 IF (contfrac_r(i,j)> epsilon) THEN 500 state_r(i,j)=is_ter 501 ELSE 502 state_r(i,j)=is_oce 503 ENDIF 504 ENDIF 505 ENDDO 506 ENDDO 507 508 509 DO j=1,nj 510 DO i=1,ni 511 IF (trip_r(i,j)>=100) THEN 512 trip_r(i,j)=norouting 513 ENDIF 514 ENDDO 515 ENDDO 516 517 CALL xios_context_initialize("orchidee_init_routing",MPI_COMM_ORCH) 518 CALL xios_orchidee_change_context("orchidee_init_routing") 519 CALL xios_set_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj, ni_glo=ni_glo, nj_glo=nj_glo) ! get routing domain dimension 520 CALL xios_set_domain_attr("routing_domain", lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! get routing domain dimension 521 CALL xios_close_context_definition() 522 CALL xios_update_calendar(1) 523 524 CALL xios_send_field("trip_r_init",trip_r) 525 526 CALL xios_send_field("trip_r",trip_r) 527 CALL xios_recv_field("trip_rp1",trip_rp1) 528 529 CALL xios_send_field("state_r",state_r) 530 CALL xios_recv_field("state_rp1",state_rp1) 531 532 ! DO j=1,nj 533 ! DO i=1,ni 534 ! IF (trip_rp1(i,j)<=8) THEN 535 ! CALL next_trip(trip_rp1, ni, nj, i, j, nexti, nextj) 536 ! IF (trip_rp1(nexti,nextj)>=100) trip_rp1(i,j)=98 ! unrouted point becomes coastal flow 537 ! ENDIF 538 ! ENDDO 539 ! ENDDO 540 541 state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj) 542 trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj) 543 544 updated=1 545 it=2 546 DO WHILE(updated>0) 547 updated=0 548 549 CALL xios_update_calendar(it) 550 551 CALL xios_send_field("trip_r", trip_r) 552 CALL xios_recv_field("trip_rp1",trip_rp1) 553 554 CALL xios_send_field("state_r",state_r) 555 CALL xios_recv_field("state_rp1",state_rp1) 556 557 state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj) 558 trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj) 559 560 ! first => riverflow and coastalflow routing to ocean must be routed to coast => to long 561 DO j=1,nj 562 DO i=1,ni 563 CALL next_trip(trip_rp1, ni, nj, i, j, nexti, nextj, ok) 564 IF (ok) THEN 565 IF (state_rp1(nexti,nextj)==is_oce .AND. (trip_rp1(nexti,nextj)==riverflow .OR. trip_rp1(nexti,nextj)==coastalflow .OR. trip_rp1(nexti,nextj)==lakeinflow)) THEN 566 trip_r(i,j) = trip_rp1(nexti,nextj) 567 updated=updated + 1 568 ENDIF 569 ENDIF 570 ENDDO 571 ENDDO 572 573 DO j=1,nj 574 DO i=1,ni 575 IF (state_rp1(i,j)==is_oce) THEN 576 IF (trip_rp1(i,j)==riverflow .OR. trip_rp1(i,j)==coastalflow .OR. trip_rp1(i,j)==lakeinflow) THEN 577 trip_r(i,j) = norouting 578 updated=updated+1 579 ENDIF 580 ENDIF 581 ENDDO 582 ENDDO 583 584 CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) 585 it = it +1 586 ENDDO 587 588 ! second prolongate riverflow and coastalflow to coast if they are in land 589 updated=1 590 DO WHILE(updated>0) 591 updated=0 592 593 CALL xios_update_calendar(it) 594 595 CALL xios_send_field("trip_r", trip_r) 596 CALL xios_recv_field("trip_rp1",trip_rp1) 597 598 CALL xios_send_field("state_r",state_r) 599 CALL xios_recv_field("state_rp1",state_rp1) 600 601 state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj) 602 trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj) 603 604 DO j=1,nj 605 DO i=1,ni 606 IF (state_rp1(i,j)==is_ter .AND. (trip_rp1(i,j)==riverflow .OR. trip_rp1(i,j)==coastalflow .OR. trip_rp1(i,j)==lakeinflow)) THEN 607 CALL next_trip(trip_extended_rp1, ni, nj, i, j, nexti, nextj, ok) 608 IF (ok) THEN 609 IF (state_rp1(nexti,nextj)==is_ter .OR. state_rp1(nexti,nextj)==is_coast) THEN 610 trip_r(i,j) = trip_extended_rp1(i,j) 611 topoind_r(i,j) = 1e-10 ! flow instantanesly for now but better to take the topind of the incomming flux 612 updated=updated+1 613 ENDIF 614 ENDIF 615 ENDIF 616 ENDDO 617 ENDDO 618 619 DO j=1,nj 620 DO i=1,ni 621 IF (state_rp1(i,j)==is_ter .OR. state_rp1(i,j)==is_coast) THEN 622 DO k=-1,1 623 DO m=-1,1 624 IF (k==0 .AND. m==0) CYCLE 625 prevj = j+k ; previ = i+m 626 CALL next_trip(trip_extended_rp1, ni, nj, previ, prevj, nexti, nextj, ok) 627 IF (ok) THEN 628 IF (nexti==i .AND. nextj==j) THEN 629 IF (state_rp1(previ,prevj)==is_ter .AND. (trip_rp1(previ,prevj)==riverflow .OR. trip_rp1(previ,prevj)==coastalflow .OR. trip_rp1(previ,prevj)==lakeinflow)) THEN 630 IF (trip_r(i,j)/=riverflow) THEN 631 trip_r(i,j)=trip_rp1(previ,prevj) 632 ELSE 633 trip_r(i,j)=riverflow 634 ENDIF 635 updated=updated+1 636 ENDIF 637 ENDIF 638 ENDIF 639 ENDDO 640 ENDDO 641 ENDIF 642 ENDDO 643 ENDDO 644 645 CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) 646 it = it +1 647 ENDDO 648 649 CALL xios_update_calendar(it) 650 651 CALL xios_send_field("trip_r", trip_r) 652 CALL xios_recv_field("trip_rp1",trip_rp1) 653 654 CALL xios_send_field("state_r",state_r) 655 CALL xios_recv_field("state_rp1",state_rp1) 656 DO j=1,nj 657 DO i=1,ni 658 IF (state_rp1(i,j)==is_coast) THEN 659 IF (trip_r(i,j)==norouting) trip_r(i,j) = coastalflow 660 ENDIF 661 ENDDO 662 ENDDO 663 664 ! it = it +1 665 ! CALL xios_update_calendar(it) 666 ! 667 ! CALL xios_send_field("trip_r", trip_r) 668 ! CALL xios_recv_field("trip_rp1",trip_rp1) 669 ! 670 ! CALL xios_send_field("state_r",state_r) 671 ! CALL xios_recv_field("state_rp1",state_rp1) 672 ! 673 ! DO j=1,nj 674 ! DO i=1,ni 675 ! IF (state_rp1(i,j)==is_oce) THEN 676 ! trip_r(i,j) = norouting 677 ! ENDIF 678 ! ENDDO 679 ! ENDDO 680 681 it = it +1 682 CALL xios_update_calendar(it) 683 684 CALL xios_send_field("trip_r", trip_r) 685 CALL xios_recv_field("trip_rp1",trip_rp1) 686 687 CALL xios_send_field("state_r",state_r) 688 CALL xios_recv_field("state_rp1",state_rp1) 689 ! check 690 updated=0 691 DO j=1,nj 692 DO i=1,ni 693 IF (state_rp1(i,j)==is_ter) THEN 694 IF (trip_r(i,j) == coastalflow .OR. trip_r(i,j) == riverflow .OR. trip_r(i,j) == norouting) THEN 695 updated=updated+1 696 trip_r(i,j) = lakeinflow 697 ENDIF 698 ELSE IF (state_rp1(i,j)==is_coast) THEN 699 IF (trip_r(i,j) == norouting) THEN 700 updated=updated+1 701 trip_r(i,j) = coastalflow 702 ENDIF 703 ENDIF 704 ENDDO 705 ENDDO 706 707 708 CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) 709 IF (updated/=0 .AND. is_mpi_root) PRINT*,"WARNING : ",updated," riverflow or coastal flow are not on the coast, ie in middle of land ==> converted into lakeinflow" 710 CALL MPI_BARRIER(MPI_COMM_ORCH, ierr) 711 712 updated=0 713 DO j=1,nj 714 DO i=1,ni 715 CALL next_trip(trip_rp1,ni,nj,i,j,nexti,nextj,ok) 716 IF (ok) THEN 717 IF (trip_rp1(nexti,nextj)==norouting) THEN 718 updated=updated+1 719 PRINT*,"point",ibegin+i,jbegin+j," (trip(i,j)=",trip_rp1(i,j),") is routed to nothing => ",ibegin+nexti,jbegin+nextj, & 720 " (trip(i,j)=",trip_rp1(nexti,nextj),")" 721 ENDIF 722 ENDIF 723 ENDDO 724 ENDDO 725 726 CALL MPI_BARRIER(MPI_COMM_ORCH, ierr) 727 CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) 728 IF (updated/=0) THEN 729 IF (is_mpi_root) PRINT*,"ERROR : ",updated," points are not correctly routed (routed to nothing)" 730 CALL xios_context_finalize 731 STOP 732 ENDIF 733 734 735 updated=0 736 DO j=1,nj 737 DO i=1,ni 738 CALL next_trip(trip_rp1,ni,nj,i,j,nexti,nextj,ok) 739 IF (ok) THEN 740 IF (state_rp1(nexti,nextj)==is_oce .AND. (trip_rp1(nexti,nextj)==riverflow .OR. trip_rp1(nexti,nextj)==coastalflow .OR. trip_rp1(nexti,nextj)==lakeinflow )) THEN 741 updated=updated+1 742 PRINT*,"point",ibegin+i,jbegin+j," (trip(i,j)=",trip_rp1(i,j),") is routed to coastalflow/riverflow/lakeinflow point that is not on land grid => ",ibegin+nexti,jbegin+nextj, & 743 " (trip(i,j)=",trip_rp1(nexti,nextj),")" 744 ENDIF 745 ENDIF 746 ENDDO 747 ENDDO 748 749 CALL MPI_BARRIER(MPI_COMM_ORCH, ierr) 750 CALL MPI_ALLREDUCE(MPI_IN_PLACE , updated, 1, MPI_INT_ORCH, MPI_SUM, MPI_COMM_ORCH, ierr) 751 IF (updated/=0) THEN 752 IF (is_mpi_root) PRINT*,"ERROR : ",updated," points are routed to coastalflow/riverflow/lakeinflow points that are not on land grid" 753 CALL xios_context_finalize 754 STOP 755 ENDIF 756 757 758 759 it = it +1 760 CALL xios_update_calendar(it) 761 762 CALL xios_send_field("trip_r", trip_r) 763 CALL xios_recv_field("trip_rp1",trip_rp1) 764 765 CALL xios_send_field("state_r",state_r) 766 CALL xios_recv_field("state_rp1",state_rp1) 767 state_r(1:ni,1:nj) = state_rp1(1:ni,1:nj) 768 trip_r(1:ni,1:nj) = trip_rp1(1:ni,1:nj) 769 770 CALL xios_send_field("trip_r_final",trip_r) 771 772 CALL xios_context_finalize 773 CALL xios_orchidee_change_context("orchidee") 774 775 CONTAINS 776 777 SUBROUTINE next_trip(trip, ni, nj, i, j, nexti, nextj, ok) 778 IMPLICIT NONE 779 REAL, INTENT(IN) :: trip(0:ni+1,0:nj+1) 780 INTEGER, INTENT(IN) :: ni 781 INTEGER, INTENT(IN) :: nj 782 INTEGER, INTENT(IN) :: i 783 INTEGER, INTENT(IN) :: j 784 INTEGER, INTENT(OUT) :: nexti 785 INTEGER, INTENT(OUT) :: nextj 786 LOGICAL, INTENT(OUT) :: ok 787 788 ok=.TRUE. 789 790 SELECT CASE(NINT(trip(i,j))) 791 CASE (1) 792 nextj=j-1 ; nexti=i ! north 793 CASE (2) 794 nextj=j-1 ; nexti=i+1 ! north-east 795 CASE (3) 796 nextj=j ; nexti=i+1 ! east 797 CASE (4) 798 nextj=j+1 ; nexti=i+1 ! south-east 799 CASE (5) 800 nextj=j+1 ; nexti=i ! south 801 CASE(6) 802 nextj=j+1 ; nexti=i-1 ! south-west 803 CASE (7) 804 nextj=j ; nexti=i-1 ! west 805 CASE (8) 806 nextj=j-1 ; nexti=i-1 ! north-west 807 CASE DEFAULT 808 ok=.FALSE. 809 END SELECT 810 811 END SUBROUTINE next_trip 812 813 END SUBROUTINE new_routing_flow_correct_riverflow 814 815 816 817 818 SUBROUTINE routing_flow_correct_riverflow(ni, nj, contfrac, coastline, trip_r, trip_extended_r, topoind_r) 819 USE xios 820 ! USE grid, ONLY : global 821 IMPLICIT NONE 822 INCLUDE "mpif.h" 823 INTEGER, INTENT(IN) :: ni ! INPUT : size i (longitude) of the local routing domain 824 INTEGER, INTENT(IN) :: nj ! INPUT : size i (latitude) of the local routing domain 825 REAL(r_std), INTENT(IN) :: contfrac(nbp_mpi) ! INPUT : continental fraction (unittless) 826 LOGICAL, INTENT(IN) :: coastline(nbp_mpi) ! INPUT/OUTPUT : coastline mask 827 REAL(r_std), INTENT(INOUT) :: trip_r(ni,nj) ! INPUT/OUTPUT : diection of flow, which will be modified by the routine 828 REAL(r_std), INTENT(INOUT) :: trip_extended_r(ni,nj) ! INPUT : direction of flow extended to ocean points 829 REAL(r_std), INTENT(INOUT) :: topoind_r(ni,nj) ! INPUT/OUTPUT : topographic index which will be modified by the routine if extended to ocean 830 831 REAL(r_std) :: contfrac_r(ni,nj) 832 REAL(r_std) :: mask_coast(nbp_mpi) 833 REAL(r_std) :: frac_coast_r(ni,nj) 834 REAL(r_std) :: trip_extended_rp1(0:ni+1,0:nj+1) 835 REAL(r_std) :: trip_rp1(0:ni+1,0:nj+1) 836 REAL(r_std) :: state_r(ni,nj) 837 REAL(r_std) :: state_rp1(0:ni+1,0:nj+1) 838 839 INTEGER, PARAMETER :: is_ter=0 840 INTEGER, PARAMETER :: is_coast=1 841 INTEGER, PARAMETER :: is_oce=2 842 INTEGER :: updated 843 INTEGER :: it, i,j,k,m, nexti,nextj,previ,prevj,ierr 844 INTEGER :: ibegin, jbegin, ni_glo, nj_glo 845 REAL(r_std) :: lonvalue_1d(ni), latvalue_1d(nj) 846 REAL(r_std),PARAMETER :: epsilon=1e-5 847 848 TYPE(xios_duration) :: ts 849 TYPE(xios_domain) :: domain_hdl 850 TYPE(xios_domaingroup) :: domain_def_hdl 851 TYPE(xios_field) :: field_hdl 852 TYPE(xios_fieldgroup) :: field_def_hdl 853 TYPE(xios_file) :: file_hdl 854 TYPE(xios_filegroup) :: file_def_hdl 855 TYPE(xios_expand_domain) :: domain_expand_hdl 856 TYPE(xios_date) :: start_date, time_origin 857 CHARACTER(LEN=20) :: calendar_type 858 LOGICAL :: true=.TRUE. 859 LOGICAL :: global=.FALSE. 422 860 423 861 CALL xios_get_domain_attr("routing_domain", ibegin=ibegin, jbegin=jbegin, ni_glo=ni_glo,nj_glo=nj_glo) ! get routing domain dimension … … 493 931 ENDDO 494 932 ENDDO 495 496 ! create a new XIOS context to solve iterative process of correction497 ! so we can use filter to update halo498 ! CALL xios_get_calendar_type(calendar_type)499 ! CALL xios_get_start_date(start_date)500 ! CALL xios_get_time_origin(time_origin)501 !502 ! CALL xios_context_initialize("orchidee_routing_correction",MPI_COMM_ORCH)503 ! CALL xios_orchidee_change_context("orchidee_routing_correction")504 ! ts.second=1505 ! calendar_type="gregorian"506 ! CALL xios_define_calendar(type=calendar_type,start_date=start_date,time_origin=time_origin, timestep=ts )507 508 ! CALL xios_get_handle("domain_definition",domain_def_hdl)509 ! CALL xios_add_child(domain_def_hdl, domain_hdl, "routing_domain")510 ! CALL xios_set_attr(domain_hdl, type="rectilinear", ibegin=ibegin, jbegin=jbegin, ni=ni, nj=nj, ni_glo=ni_glo, nj_glo=nj_glo) ! set routing domain dimension511 ! CALL xios_set_attr(domain_hdl, lonvalue_1d=lonvalue_1d, latvalue_1d=latvalue_1d) ! set routing domain dimension512 ! CALL xios_add_child(domain_def_hdl, domain_hdl, "routing_domain_expand")513 ! CALL xios_set_attr(domain_hdl, domain_ref="routing_domain")514 ! CALL xios_add_child(domain_hdl,domain_expand_hdl)515 ! CALL xios_set_attr(domain_expand_hdl, type="edge", i_periodic=.TRUE., j_periodic=.TRUE.)516 !517 ! CALL xios_get_handle("field_definition",field_def_hdl)518 ! CALL xios_add_child(field_def_hdl,field_hdl,"trip_r_init")519 ! CALL xios_set_attr(field_hdl,domain_ref="routing_domain")520 ! CALL xios_add_child(field_def_hdl,field_hdl,"trip_r_final")521 ! CALL xios_set_attr(field_hdl,domain_ref="routing_domain")522 523 ! CALL xios_add_child(field_def_hdl,field_hdl,"trip_r")524 ! CALL xios_set_attr(field_hdl,domain_ref="routing_domain")525 ! CALL xios_add_child(field_def_hdl,field_hdl,"trip_rp1")526 ! CALL xios_set_attr(field_hdl, field_ref="trip_r", domain_ref="routing_domain_expand", read_access=.TRUE.)527 528 ! CALL xios_add_child(field_def_hdl,field_hdl,"state_r")529 ! CALL xios_set_attr(field_hdl,domain_ref="routing_domain")530 ! CALL xios_add_child(field_def_hdl,field_hdl,"state_rp1")531 ! CALL xios_set_attr(field_hdl, field_ref="state_r", domain_ref="routing_domain_expand", read_access=.TRUE.)532 !533 ! CALL xios_get_handle("file_definition",file_def_hdl)534 ! CALL xios_add_child(file_def_hdl,file_hdl,"routing_correction")535 ! CALL xios_set_attr(file_hdl,type="one_file", output_freq=ts, sync_freq=ts, enabled=.TRUE.)536 ! CALL xios_add_child(file_hdl, field_hdl)537 ! CALL xios_set_attr(field_hdl, field_ref="trip_r", operation="instant")538 ! CALL xios_add_child(file_hdl, field_hdl)539 ! CALL xios_set_attr(field_hdl, field_ref="state_r", operation="instant")540 ! CALL xios_add_child(file_hdl, field_hdl)541 ! CALL xios_set_attr(field_hdl, field_ref="trip_r_init", operation="once")542 ! CALL xios_add_child(file_hdl, field_hdl)543 ! CALL xios_set_attr(field_hdl, field_ref="state_r_init", operation="once")544 545 933 546 934 CALL xios_context_initialize("orchidee_init_routing",MPI_COMM_ORCH) … … 831 1219 LOGICAL :: file_exists 832 1220 REAL(r_std) :: epsilon = 1e-5 833 1221 CHARACTER(LEN=255) :: routing_file_type 834 1222 !_ ================================================================================================================================ 835 1223 ! … … 879 1267 !Config Units = [days] 880 1268 1269 CALL getin_p('SLOW_TCST', slow_tcst) 1270 1271 routing_file_type = "standard" 1272 CALL getin_p("routing_file_type", routing_file_type) 1273 1274 IF (TRIM(routing_file_type)=="standard") THEN 1275 topoind_factor=1000 1276 ELSE IF (TRIM(routing_file_type)=="merit") THEN 1277 topoind_factor=1 1278 ELSE 1279 CALL ipslerr_p(3,'routing_flow_init_local', & 1280 'Error in call routing_flow_init_local','getin routing_file_type : bad value, must be <standard> or <merit>','') 1281 ENDIF 1282 881 1283 split_routing=1 882 1284 CALL getin_p("SPLIT_ROUTING",split_routing) … … 892 1294 893 1295 IF (is_omp_root) THEN 894 895 WHERE(contfrac_mpi<=epsilon) contfrac_mpi=0 896 WHERE(contfrac_mpi>=epsilon) contfrac_mpi=1 897 898 diag=0 899 WHERE(coastline_mpi) diag=1 900 CALL xios_send_field("is_coastline", diag) 901 1296 WHERE(contfrac_mpi <= epsilon) contfrac_mpi=0 1297 WHERE(contfrac_mpi >= 1-epsilon) contfrac_mpi=1 1298 902 1299 CALL xios_get_domain_attr("routing_domain", ni=ni, nj=nj) ! get routing domain dimension 903 1300 … … 933 1330 934 1331 ALLOCATE(trip_extended_r(ni,nj)) 1332 1333 ! correction on coastline in case of LAM, the cells near the frontier are considered as coastline 1334 ! => interpolate model grid to routing grid and reinterpolate to model grid. Fractionnal cells will be considered as coastline too 1335 ! diag=1 1336 ! CALL xios_send_field("one", diag) 1337 ! CALL xios_recv_field("tmp1_coastline", diag_r) 1338 ! CALL xios_send_field("tmp2_coastline", diag_r) 1339 ! CALL xios_recv_field("lam_coastline", diag) 1340 ! WHERE(diag<1-epsilon) is_coastline=.TRUE. 1341 1342 diag=0 1343 WHERE(is_coastline) diag=1 1344 CALL xios_send_field("is_coastline", diag) 935 1345 936 1346 is_lakeinflow_r(:) = .FALSE. … … 960 1370 CALL xios_recv_field("basins_extended_r",basins_extended_r) ! recv basins index from file 961 1371 962 CALL routing_flow_correct_riverflow(ni, nj, contfrac_mpi, coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r) 1372 ! CALL routing_flow_correct_riverflow(ni, nj, contfrac_mpi, coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r) 1373 CALL new_routing_flow_correct_riverflow(ni, nj, contfrac_mpi, is_coastline, trip_rp1(1:ni,1:nj), trip_extended_r, topoind_r) 963 1374 964 1375 CALL xios_send_field("trip_update_r",trip_rp1(1:ni,1:nj)) ! send to xios trip array to update for halo … … 1041 1452 CALL xios_send_field("routing_mask_r", diag_r) 1042 1453 CALL xios_recv_field("frac_routing", diag) 1454 WHERE (diag < epsilon) diag=0 1043 1455 routing_weight_in (:) = 0 1044 1456 unrouted_weight (:) = 0 … … 1049 1461 ! compute the number of coast cells to distribute lakeinflow onto 1050 1462 diag=0 1051 WHERE ( coastline_mpi) diag=11463 WHERE (is_coastline) diag=1 1052 1464 nb_coast_points=SUM(diag) 1053 1465 CALL reduce_sum_mpi(nb_coast_points, total_coast_points) … … 1067 1479 1068 1480 diag(:) = 0 1069 WHERE ( coastline_mpi) diag=11481 WHERE (is_coastline) diag=1 1070 1482 CALL xios_send_field("mask_coastal",diag) 1071 1483 CALL xios_recv_field("frac_coastal_r",diag_r) 1072 WHERE (diag_r > 100 .) ! missing_value1484 WHERE (diag_r > 100) ! missing_value 1073 1485 diag_r=0 1074 1486 ELSE WHERE (diag_r < epsilon) … … 1080 1492 1081 1493 diag(:) = 0 1082 WHERE (.NOT. coastline_mpi) diag=11494 WHERE (.NOT. is_coastline) diag=1 1083 1495 CALL xios_send_field("mask_lake",diag) 1084 1496 CALL xios_recv_field("frac_lake_r",diag_r) 1085 WHERE (diag_r > 100 .) ! missing_value1497 WHERE (diag_r > 100) ! missing_value 1086 1498 diag_r=0. 1087 1499 ELSEWHERE (diag_r < epsilon) … … 1105 1517 weight_coast_to_lake_r(ij) = 1./frac_lake_r(ij) 1106 1518 weight_coast_to_coast_r(ij) = 0 1107 ELSE IF (frac_coast_r(ij)>0 1108 weight_coast_to_coast_r(ij) = 1./(frac_coast_r(ij) + frac_lake_r(ij))1519 ELSE IF (frac_coast_r(ij)>0) THEN 1520 weight_coast_to_coast_r(ij) = 1./(frac_coast_r(ij)) 1109 1521 weight_coast_to_lake_r(ij) = 0 1522 ELSE 1110 1523 ENDIF 1111 1524 ELSE IF (is_lakeinflow_r(ij)) THEN … … 1117 1530 weight_lake_to_lake_r(ij) = 0 1118 1531 ELSE IF (frac_lake_r(ij)>0 ) THEN 1119 weight_lake_to_lake_r(ij) = 1./(frac_ coast_r(ij) + frac_lake_r(ij))1532 weight_lake_to_lake_r(ij) = 1./(frac_lake_r(ij)) 1120 1533 weight_lake_to_coast_r(ij) = 0 1121 1534 ENDIF … … 1432 1845 WHERE(.NOT. routing_mask_r) drainage_r = 0 1433 1846 1434 CALL MPI_ALLREDUCE(sum(runoff* routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)1847 CALL MPI_ALLREDUCE(sum(runoff*(routing_weight-unrouted_weight)),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 1435 1848 CALL MPI_ALLREDUCE(sum(runoff_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 1436 1849 IF (sum_water_after/=0) THEN … … 1441 1854 ENDIF 1442 1855 1443 CALL MPI_ALLREDUCE(sum(drainage* routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)1856 CALL MPI_ALLREDUCE(sum(drainage*(routing_weight-unrouted_weight)),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 1444 1857 CALL MPI_ALLREDUCE(sum(drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 1445 1858 IF (sum_water_after/=0) THEN … … 1451 1864 1452 1865 1453 CALL MPI_ALLREDUCE(sum( runoff*routing_weight)+sum(drainage*routing_weight),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)1454 CALL MPI_ALLREDUCE(sum(runoff_r ) +sum(drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr)1866 CALL MPI_ALLREDUCE(sum((runoff+drainage)*(routing_weight-unrouted_weight)),sum_water_before,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 1867 CALL MPI_ALLREDUCE(sum(runoff_r+drainage_r),sum_water_after,1,MPI_REAL_ORCH,MPI_SUM,MPI_COMM_ORCH,ierr) 1455 1868 1456 1869 IF (is_mpi_root) PRINT *,"Loose water by interpolation ; before : ", sum_water_before," ; after : ",sum_water_after, & … … 1470 1883 IF ( routing_mask_r(ig) ) THEN 1471 1884 1472 flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig) /1000.)*fast_tcst*one_day/(dt_routing/split_routing)),&1885 flow = MIN(fast_reservoir_r(ig)/((topoind_r(ig)*topoind_factor/1000.)*fast_tcst*one_day/(dt_routing/split_routing)),& 1473 1886 & fast_reservoir_r(ig)-min_sechiba) 1474 1887 fast_flow_r(ig) = MAX(flow, zero) 1475 1888 ! 1476 flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig) /1000.)*slow_tcst*one_day/(dt_routing/split_routing)),&1889 flow = MIN(slow_reservoir_r(ig)/((topoind_r(ig)*topoind_factor/1000.)*slow_tcst*one_day/(dt_routing/split_routing)),& 1477 1890 & slow_reservoir_r(ig)-min_sechiba) 1478 1891 slow_flow_r(ig) = MAX(flow, zero) 1479 1892 ! 1480 flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig) /1000.)*stream_tcst*one_day/(dt_routing/split_routing)),&1893 flow = MIN(stream_reservoir_r(ig)/((topoind_r(ig)*topoind_factor/1000.)*stream_tcst*one_day/(dt_routing/split_routing)),& 1481 1894 & stream_reservoir_r(ig)-min_sechiba) 1482 1895 stream_flow_r(ig) = MAX(flow, zero) … … 1641 2054 ENDIF 1642 2055 1643 coastalflow=coastalflow+flow_coast 2056 coastalflow=coastalflow+flow_coast+(runoff+drainage)*unrouted_weight 1644 2057 lakeinflow=lakeinflow+flow_lake 1645 2058 1646 2059 1647 WHERE(is_coastline) coastalflow = coastalflow + runoff_in + drainage_in1648 WHERE(.NOT.is_coastline) lakeinflow = lakeinflow + runoff_in + drainage_in2060 ! WHERE(is_coastline) coastalflow = coastalflow + runoff_in + drainage_in 2061 ! WHERE(.NOT.is_coastline) lakeinflow = lakeinflow + runoff_in + drainage_in 1649 2062 1650 2063 sum_water_after = sum(coastalflow(:) + riverflow(:) + lakeinflow(:))+sum(fast_reservoir_r(:)+slow_reservoir_r(:)+stream_reservoir_r(:))
Note: See TracChangeset
for help on using the changeset viewer.
