Changeset 4816
- Timestamp:
- Feb 13, 2024, 6:17:52 PM (3 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lmdz_wake.F90
r4814 r4816 227 227 228 228 ! Variables liees au calcul de hw 229 REAL, DIMENSION (klon) :: ptop _provis, ptop, ptop_new229 REAL, DIMENSION (klon) :: ptop 230 230 REAL, DIMENSION (klon) :: sum_dth 231 231 REAL, DIMENSION (klon) :: dthmin … … 311 311 312 312 LOGICAL :: first_call=.true. 313 313 314 314 315 ! ------------------------------------------------------------------------- … … 697 698 ! Determine Wake top pressure (Ptop) from buoyancy integral 698 699 ! -------------------------------------------------------- 699 700 ! -1/ Pressure of the level where dth becomes less than delta_t_min. 701 702 DO i = 1, klon 703 ptop_provis(i) = ph(i, 1) 704 END DO 705 DO k = 2, klev 706 DO i = 1, klon 707 708 ! IM v3JYG; ptop_provis(i).LT. ph(i,1) 709 710 ! LJYF : mettre cette ligne à la place de la suivante IF (dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. & 711 IF (dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min .AND. & 712 ptop_provis(i)==ph(i,1)) THEN 713 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)- & 714 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 715 END IF 716 END DO 717 END DO 718 719 ! -2/ dth integral 720 721 DO i = 1, klon 722 sum_dth(i) = 0. 723 dthmin(i) = -delta_t_min 724 z(i) = 0. 725 END DO 726 727 DO k = 1, klev 728 DO i = 1, klon 729 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG) 730 IF (dz(i)>0) THEN 731 z(i) = z(i) + dz(i) 732 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 733 dthmin(i) = amin1(dthmin(i), dth(i,k)) 734 END IF 735 END DO 736 END DO 737 738 ! -3/ height of triangle with area= sum_dth and base = dthmin 739 740 DO i = 1, klon 741 hw0(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) 742 hw0(i) = amax1(hwmin, hw0(i)) 743 END DO 744 745 ! -4/ now, get Ptop 746 747 DO i = 1, klon 748 z(i) = 0. 749 ptop(i) = ph(i, 1) 750 END DO 751 752 DO k = 1, klev 753 DO i = 1, klon 754 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw0(i)-z(i)) 755 IF (dz(i)>0) THEN 756 z(i) = z(i) + dz(i) 757 ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i) 758 END IF 759 END DO 760 END DO 761 762 IF (prt_level>=10) THEN 763 PRINT *, 'wake-2, ptop_provis(igout), ptop(igout) ', ptop_provis(igout), ptop(igout) 764 ENDIF 765 766 767 ! -5/ Determination de ktop et kupper 768 769 DO k = klev, 1, -1 770 DO i = 1, klon 771 IF (ph(i,k+1)<ptop(i)) ktop(i) = k 772 END DO 773 END DO 774 !print*, 'ptop, pupper, ktop, kupper', ptop, pupper, ktop, kupper 775 776 777 778 ! -6/ Correct ktop and ptop 779 780 DO i = 1, klon 781 ptop_new(i) = ptop(i) 782 END DO 783 DO k = klev, 2, -1 784 DO i = 1, klon 785 IF (k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & 786 ! LJYF changer par cette ligne dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 787 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 788 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1)-(dth(i, & 789 k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 790 END IF 791 END DO 792 END DO 793 794 DO i = 1, klon 795 ptop(i) = ptop_new(i) 796 END DO 797 798 DO k = klev, 1, -1 799 DO i = 1, klon 800 IF (ph(i,k+1)<ptop(i)) ktop(i) = k 801 END DO 802 END DO 803 804 CALL pkupper (klon, klev, ptop, ph, pupper, kupper) 805 806 IF (prt_level>=10) THEN 807 PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) 700 Do i=1, klon 701 wk_adv(i) = .True. 702 Enddo 703 Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, & 704 dth, hw0, rho, & 705 ktop, wk_adv) 706 707 IF (prt_level>=10) THEN 708 PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) 808 709 ENDIF 809 710 … … 848 749 dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG) 849 750 IF (dz(i)>0) THEN 751 ! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau 850 752 z(i) = z(i) + dz(i) 851 753 sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i) … … 988 890 ! ------------------------------------------------------------------------ 989 891 ! 892 wk_adv = .True. 990 893 DO isubstep = 1, nsub 991 894 ! … … 993 896 ! wk_adv is the logical flag enabling wake evolution in the time advance 994 897 ! loop 898 995 899 DO i = 1, klon 996 900 wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1. … … 1826 1730 1827 1731 1828 ! Determine Ptop from buoyancy integral 1829 ! --------------------------------------- 1830 1831 ! - 1/ Pressure of the level where dth changes sign. 1832 1833 DO i = 1, klon 1834 IF (wk_adv(i)) THEN 1835 ptop_provis(i) = ph(i, 1) 1836 END IF 1837 END DO 1838 1839 DO k = 2, klev 1840 DO i = 1, klon 1841 IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & 1842 ! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 1843 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 1844 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 1845 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 1846 END IF 1847 END DO 1848 END DO 1849 1850 ! - 2/ dth integral 1851 1852 DO i = 1, klon 1853 IF (wk_adv(i)) THEN !!! nrlmd 1854 sum_dth(i) = 0. 1855 dthmin(i) = -delta_t_min 1856 z(i) = 0. 1857 END IF 1858 END DO 1859 1860 DO k = 1, klev 1861 DO i = 1, klon 1862 IF (wk_adv(i)) THEN 1863 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG) 1864 IF (dz(i)>0) THEN 1865 z(i) = z(i) + dz(i) 1866 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 1867 dthmin(i) = amin1(dthmin(i), dth(i,k)) 1868 END IF 1869 END IF 1870 END DO 1871 END DO 1872 1873 ! - 3/ height of triangle with area= sum_dth and base = dthmin 1874 1875 DO i = 1, klon 1876 IF (wk_adv(i)) THEN 1877 hw(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) 1878 hw(i) = amax1(hwmin, hw(i)) 1879 END IF 1880 END DO 1881 1882 ! - 4/ now, get Ptop 1883 1884 DO i = 1, klon 1885 IF (wk_adv(i)) THEN !!! nrlmd 1886 ktop(i) = 0 1887 z(i) = 0. 1888 END IF 1889 END DO 1890 1891 DO k = 1, klev 1892 DO i = 1, klon 1893 IF (wk_adv(i)) THEN 1894 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw(i)-z(i)) 1895 IF (dz(i)>0) THEN 1896 z(i) = z(i) + dz(i) 1897 ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i) 1898 ktop(i) = k 1899 END IF 1900 END IF 1901 END DO 1902 END DO 1903 1904 ! 4.5/Correct ktop and ptop 1905 1906 DO i = 1, klon 1907 IF (wk_adv(i)) THEN 1908 ptop_new(i) = ptop(i) 1909 END IF 1910 END DO 1911 1912 DO k = klev, 2, -1 1913 DO i = 1, klon 1914 ! IM v3JYG; IF (k .GE. ktop(i) 1915 IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & 1916 ! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 1917 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 1918 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 1919 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 1920 END IF 1921 END DO 1922 END DO 1923 1924 1925 DO i = 1, klon 1926 IF (wk_adv(i)) THEN 1927 ptop(i) = ptop_new(i) 1928 END IF 1929 END DO 1930 1931 DO k = klev, 1, -1 1932 DO i = 1, klon 1933 IF (wk_adv(i)) THEN !!! nrlmd 1934 IF (ph(i,k+1)<ptop(i)) ktop(i) = k 1935 END IF 1936 END DO 1937 END DO 1938 1939 CALL pkupper (klon, klev, ptop, ph, pupper, kupper) 1732 1733 Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, & 1734 dth, hw, rho, & 1735 ktop, wk_adv) 1940 1736 1941 1737 … … 2651 2447 2652 2448 2653 SUBROUTINE pkupper (klon, klev, ptop, ph, pupper, kupper) 2449 SUBROUTINE pkupper (klon, klev, ptop, ph, p, pupper, kupper, & 2450 dth, hw_, rho, & 2451 ktop, wk_adv) 2654 2452 2655 2453 USE lmdz_wake_ini , ONLY : wk_pupper 2454 USE lmdz_wake_ini , ONLY : RG 2455 USE lmdz_wake_ini , ONLY : hwmin 2656 2456 IMPLICIT NONE 2657 2457 2658 2458 INTEGER, INTENT(IN) :: klon,klev 2659 REAL, INTENT(IN), DIMENSION (klon,klev+1) :: ph 2660 REAL, INTENT(IN), DIMENSION (klon) :: ptop 2661 REAL, INTENT(OUT), DIMENSION (klon) :: pupper 2662 INTEGER, INTENT(OUT), DIMENSION (klon) :: kupper 2459 REAL, INTENT(IN), DIMENSION (klon,klev+1) :: ph, p 2460 REAL, INTENT(IN), DIMENSION (klon,klev+1) :: rho 2461 LOGICAL, INTENT(IN), DIMENSION (klon) :: wk_adv 2462 REAL, INTENT(IN), DIMENSION (klon,klev+1) :: dth 2463 2464 REAL, INTENT(OUT), DIMENSION (klon) :: hw_ 2465 REAL, INTENT(OUT), DIMENSION (klon) :: ptop 2466 INTEGER, INTENT(OUT), DIMENSION (klon) :: Ktop 2467 REAL, INTENT(OUT), DIMENSION (klon) :: pupper 2468 INTEGER, INTENT(OUT), DIMENSION (klon) :: kupper 2663 2469 INTEGER :: i,k 2664 2470 REAL :: delta_t_min 2471 2472 2473 REAL, DIMENSION (klon) :: dthmin 2474 REAL, DIMENSION (klon) :: ptop_provis,ptop_new 2475 REAL, DIMENSION (klon) :: z, dz 2476 REAL, DIMENSION (klon) :: sum_dth 2477 2478 ! LJYF : a priori z, dz sum_dth sont aussi des variables internes 2479 ! Les eliminer apres verification convergence numerique 2480 2481 delta_t_min = 0.2 2482 2483 ! Determine Ptop from buoyancy integral 2484 ! --------------------------------------- 2485 2486 ! - 1/ Pressure of the level where dth changes sign. 2487 print*,'WAKE LJYF' 2488 2489 DO i = 1, klon 2490 IF (wk_adv(i)) THEN 2491 ptop_provis(i) = ph(i, 1) 2492 END IF 2493 END DO 2494 2495 DO k = 2, klev 2496 DO i = 1, klon 2497 IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & 2498 ! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 2499 dth(i,k)>-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 2500 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 2501 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 2502 END IF 2503 END DO 2504 END DO 2505 2506 ! - 2/ dth integral 2507 2508 DO i = 1, klon 2509 IF (wk_adv(i)) THEN !!! nrlmd 2510 sum_dth(i) = 0. 2511 dthmin(i) = -delta_t_min 2512 z(i) = 0. 2513 END IF 2514 END DO 2515 2516 DO k = 1, klev 2517 DO i = 1, klon 2518 IF (wk_adv(i)) THEN 2519 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG) 2520 IF (dz(i)>0) THEN 2521 z(i) = z(i) + dz(i) 2522 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 2523 dthmin(i) = amin1(dthmin(i), dth(i,k)) 2524 END IF 2525 END IF 2526 END DO 2527 END DO 2528 2529 ! - 3/ height of triangle with area= sum_dth and base = dthmin 2530 2531 DO i = 1, klon 2532 IF (wk_adv(i)) THEN 2533 hw_(i) = 2.*sum_dth(i)/amin1(dthmin(i), -0.5) 2534 hw_(i) = amax1(hwmin, hw_(i)) 2535 END IF 2536 END DO 2537 2538 ! - 4/ now, get Ptop 2539 2540 DO i = 1, klon 2541 IF (wk_adv(i)) THEN !!! nrlmd 2542 ktop(i) = 0 2543 z(i) = 0. 2544 END IF 2545 END DO 2546 2547 DO k = 1, klev 2548 DO i = 1, klon 2549 IF (wk_adv(i)) THEN 2550 dz(i) = amin1(-(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG), hw_(i)-z(i)) 2551 IF (dz(i)>0) THEN 2552 z(i) = z(i) + dz(i) 2553 ptop(i) = ph(i, k) - rho(i, k)*RG*dz(i) 2554 ktop(i) = k 2555 END IF 2556 END IF 2557 END DO 2558 END DO 2559 2560 ! 4.5/Correct ktop and ptop 2561 2562 DO i = 1, klon 2563 IF (wk_adv(i)) THEN 2564 ptop_new(i) = ptop(i) 2565 END IF 2566 END DO 2567 2568 DO k = klev, 2, -1 2569 DO i = 1, klon 2570 ! IM v3JYG; IF (k .GE. ktop(i) 2571 IF (wk_adv(i) .AND. k<=ktop(i) .AND. ptop_new(i)==ptop(i) .AND. & 2572 ! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 2573 dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 2574 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 2575 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 2576 END IF 2577 END DO 2578 END DO 2579 2580 2581 DO i = 1, klon 2582 IF (wk_adv(i)) THEN 2583 ptop(i) = ptop_new(i) 2584 END IF 2585 END DO 2586 2587 DO k = klev, 1, -1 2588 DO i = 1, klon 2589 IF (wk_adv(i)) THEN !!! nrlmd 2590 IF (ph(i,k+1)<ptop(i)) ktop(i) = k 2591 END IF 2592 END DO 2593 END DO 2594 2595 ! IF (prt_level>=10) THEN 2596 ! PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) 2597 ! ENDIF 2665 2598 2666 2599 kupper = 0
Note: See TracChangeset
for help on using the changeset viewer.