Changeset 4816


Ignore:
Timestamp:
Feb 13, 2024, 6:17:52 PM (3 months ago)
Author:
fhourdin
Message:

Nettoyage wake (Lamine, Jean-Yves, Frederic)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/lmdz_wake.F90

    r4814 r4816  
    227227
    228228  ! Variables liees au calcul de hw
    229   REAL, DIMENSION (klon)                                :: ptop_provis, ptop, ptop_new
     229  REAL, DIMENSION (klon)                                :: ptop
    230230  REAL, DIMENSION (klon)                                :: sum_dth
    231231  REAL, DIMENSION (klon)                                :: dthmin
    … …  
    311311
    312312  LOGICAL                                               :: first_call=.true.
     313
    313314
    314315  ! -------------------------------------------------------------------------
    … …  
    697698  ! Determine Wake top pressure (Ptop) from buoyancy integral
    698699  ! --------------------------------------------------------
    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)
    808709  ENDIF
    809710
    … …  
    848749      dz(i) = -(amax1(ph(i,k+1),ptop(i))-ph(i,k))/(rho(i,k)*RG)
    849750      IF (dz(i)>0) THEN
     751              ! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau
    850752        z(i) = z(i) + dz(i)
    851753        sum_thu(i) = sum_thu(i) + thu(i, k)*dz(i)
    … …  
    988890  ! ------------------------------------------------------------------------
    989891  !
     892  wk_adv = .True.
    990893  DO isubstep = 1, nsub
    991894  !
    … …  
    993896    ! wk_adv is the logical flag enabling wake evolution in the time advance
    994897    ! loop
     898   
    995899    DO i = 1, klon
    996900      wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1.
    … …  
    18261730
    18271731
    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)
    19401736 
    19411737
    … …  
    26512447
    26522448
    2653 SUBROUTINE pkupper (klon, klev, ptop, ph, pupper, kupper)
     2449SUBROUTINE pkupper (klon, klev, ptop, ph, p, pupper, kupper, &
     2450                    dth, hw_, rho,  &
     2451                    ktop, wk_adv)
    26542452
    26552453USE lmdz_wake_ini , ONLY : wk_pupper
     2454USE lmdz_wake_ini , ONLY : RG
     2455USE lmdz_wake_ini , ONLY : hwmin
    26562456IMPLICIT NONE
    26572457
    26582458INTEGER,  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
     2459REAL,     INTENT(IN),   DIMENSION (klon,klev+1)   :: ph, p
     2460REAL,     INTENT(IN),   DIMENSION (klon,klev+1)   :: rho
     2461LOGICAL,  INTENT(IN),   DIMENSION (klon)          :: wk_adv
     2462REAL,     INTENT(IN),   DIMENSION (klon,klev+1)   :: dth
     2463
     2464REAL,     INTENT(OUT),   DIMENSION (klon)         :: hw_
     2465REAL,     INTENT(OUT),   DIMENSION (klon)         :: ptop
     2466INTEGER,  INTENT(OUT),   DIMENSION (klon)         :: Ktop
     2467REAL,     INTENT(OUT),   DIMENSION (klon)         :: pupper
     2468INTEGER,  INTENT(OUT),   DIMENSION (klon)         :: kupper
    26632469INTEGER                                           :: i,k
    2664 
     2470REAL                                              :: delta_t_min
     2471
     2472
     2473REAL,     DIMENSION (klon)         :: dthmin
     2474REAL,     DIMENSION (klon)         :: ptop_provis,ptop_new
     2475REAL,     DIMENSION (klon)         :: z, dz
     2476REAL,     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
    26652598
    26662599 kupper = 0
Note: See TracChangeset for help on using the changeset viewer.