Changeset 4834
- Timestamp:
- Feb 29, 2024, 2:25:22 AM (2 months ago)
- File:
-
- 1 edited
-
LMDZ6/trunk/libf/phylmd/lmdz_wake.F90 (modified) (55 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lmdz_wake.F90
r4825 r4834 6 6 7 7 SUBROUTINE wake(klon,klev,znatsurf, p, ph, pi, dtime, & 8 t env0, qe0, omgb, &8 tb0, qb0, omgb, & 9 9 dtdwn, dqdwn, amdwn, amup, dta, dqa, wgen, & 10 10 sigd_con, Cin, & 11 11 deltatw, deltaqw, sigmaw, asigmaw, wdens, awdens, & ! state variables 12 12 dth, hw, wape, fip, gfl, & 13 dtls, dqls, ktopw, omgbdth, dp_omgb, t u, qu, &13 dtls, dqls, ktopw, omgbdth, dp_omgb, tx, qx, & 14 14 dtke, dqke, omg, dp_deltomg, wkspread, cstar, & 15 15 d_deltat_gw, & ! tendencies … … 80 80 81 81 ! aire : aire de la maille 82 ! t env0 : temperature dans l'environnement(K)83 ! q e0 : humidite dans l'environnement(kg/kg)82 ! tb0 : horizontal average of temperature (K) 83 ! qb0 : horizontal average of humidity (kg/kg) 84 84 ! omgb : vitesse verticale moyenne sur la maille (Pa/s) 85 85 ! dtdwn: source de chaleur due aux descentes (K/s) … … 101 101 ! Variables internes : 102 102 103 ! rhow : masse volumique de la poche froide 104 ! rho : environment density at P levels 105 ! rhoh : environment density at Ph levels 106 ! tenv : environment temperature | may change within 107 ! qe : environment humidity | sub-time-stepping 108 ! the : environment potential temperature 109 ! thu : potential temperature in undisturbed area 110 ! tu : temperature in undisturbed area 111 ! qu : humidity in undisturbed area 103 ! rho : mean density at P levels 104 ! rhoh : mean density at Ph levels 105 ! tb : mean temperature | may change within 106 ! qb : mean humidity | sub-time-stepping 107 ! thb : mean potential temperature 108 ! thx : potential temperature in (x) area 109 ! tx : temperature in (x) area 110 ! qx : humidity in (x) area 112 111 ! dp_omgb: vertical gradient og LS omega 113 112 ! omgbw : wake average vertical omega … … 115 114 ! omgbdq : flux of Delta_q transported by LS omega 116 115 ! dth : potential temperature diff. wake-undist. 117 ! th1 : first pot. temp. for vertical advection (=th u)116 ! th1 : first pot. temp. for vertical advection (=thx) 118 117 ! th2 : second pot. temp. for vertical advection (=thw) 119 118 ! q1 : first humidity for vertical advection 120 119 ! q2 : second humidity for vertical advection 121 ! d_deltatw : terme de redistribution pour deltatw122 ! d_deltaqw : terme de redistribution pour deltaqw123 ! deltatw0 : deltatw initial124 ! deltaqw0 : deltaqw initial120 ! d_deltatw : redistribution term for deltatw 121 ! d_deltaqw : redistribution term for deltaqw 122 ! deltatw0 : initial deltatw 123 ! deltaqw0 : initial deltaqw 125 124 ! hw0 : wake top hight (defined as the altitude at which deltatw=0) 126 125 ! amflux : horizontal mass flux through wake boundary 127 126 ! wdens_ref: initial number of wakes per unit area (3D) or per 128 ! unit length (2D), at the beginning of each time step127 ! unit length (2D), at the beginning of each time step 129 128 ! Tgw : 1 sur la periode de onde de gravite 130 129 ! Cgw : vitesse de propagation de onde de gravite 131 ! LL : distance entre 2 poches130 ! LL : distance between 2 wakes 132 131 133 132 ! ------------------------------------------------------------------------- … … 145 144 REAL, DIMENSION (klon, klev), INTENT(IN) :: omgb 146 145 REAL, INTENT(IN) :: dtime 147 REAL, DIMENSION (klon, klev), INTENT(IN) :: t env0, qe0146 REAL, DIMENSION (klon, klev), INTENT(IN) :: tb0, qb0 148 147 REAL, DIMENSION (klon, klev), INTENT(IN) :: dtdwn, dqdwn 149 148 REAL, DIMENSION (klon, klev), INTENT(IN) :: amdwn, amup … … 166 165 167 166 REAL, DIMENSION (klon, klev), INTENT(OUT) :: dth 168 REAL, DIMENSION (klon, klev), INTENT(OUT) :: t u, qu167 REAL, DIMENSION (klon, klev), INTENT(OUT) :: tx, qx 169 168 REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtls, dqls 170 169 REAL, DIMENSION (klon, klev), INTENT(OUT) :: dtke, dqke … … 195 194 REAL, DIMENSION (klon, klev) :: deltatw0 196 195 REAL, DIMENSION (klon, klev) :: deltaqw0 197 REAL, DIMENSION (klon, klev) :: t env, qe196 REAL, DIMENSION (klon, klev) :: tb, qb 198 197 199 198 ! Variables liees a la dynamique de population 1 … … 239 238 ! Sub-timestep tendencies and related variables 240 239 REAL, DIMENSION (klon, klev) :: d_deltatw, d_deltaqw 241 REAL, DIMENSION (klon, klev) :: d_t env, d_qe240 REAL, DIMENSION (klon, klev) :: d_tb, d_qb 242 241 REAL, DIMENSION (klon) :: d_wdens, d_awdens, d_sigmaw, d_asigmaw 243 242 REAL, DIMENSION (klon) :: d_sig_gen, d_sig_death, d_sig_col, d_sig_spread, d_sig_bnd … … 261 260 REAL :: d_wdens_targ 262 261 263 REAL, DIMENSION (klon) :: sum_th u, sum_tu, sum_qu, sum_thvu264 REAL, DIMENSION (klon) :: sum_dq , sum_rho262 REAL, DIMENSION (klon) :: sum_thx, sum_tx, sum_qx, sum_thvx 263 REAL, DIMENSION (klon) :: sum_dq 265 264 REAL, DIMENSION (klon) :: sum_dtdwn, sum_dqdwn 266 REAL, DIMENSION (klon) :: av_th u, av_tu, av_qu, av_thvu267 REAL, DIMENSION (klon) :: av_dth, av_dq , av_rho265 REAL, DIMENSION (klon) :: av_thx, av_tx, av_qx, av_thvx 266 REAL, DIMENSION (klon) :: av_dth, av_dq 268 267 REAL, DIMENSION (klon) :: av_dtdwn, av_dqdwn 269 268 270 REAL, DIMENSION (klon, klev) :: rho , rhow269 REAL, DIMENSION (klon, klev) :: rho 271 270 REAL, DIMENSION (klon, klev+1) :: rhoh 272 REAL, DIMENSION (klon, klev) :: rhow_moyen273 271 REAL, DIMENSION (klon, klev) :: zh 274 272 REAL, DIMENSION (klon, klev+1) :: zhh 275 REAL, DIMENSION (klon, klev) :: epaisseur1, epaisseur2 276 277 REAL, DIMENSION (klon, klev) :: the, thu 273 274 REAL, DIMENSION (klon, klev) :: thb, thx 278 275 279 276 REAL, DIMENSION (klon, klev) :: omgbw … … 312 309 LOGICAL :: first_call=.true. 313 310 311 312 !!-- variables liees au nouveau calcul de ptop et hw 313 REAL, DIMENSION (klon, klev) :: int_dth 314 REAL, DIMENSION (klon, klev) :: zzz, dzzz 315 REAL :: epsil 316 REAL, DIMENSION (klon) :: ptop1 317 INTEGER, DIMENSION (klon) :: ktop1 318 REAL, DIMENSION (klon) :: omega 319 REAL, DIMENSION (klon) :: h_zzz 320 321 ! print*,'WAKE LJYF' 314 322 315 323 ! ------------------------------------------------------------------------- … … 388 396 389 397 delta_t_min = 0.2 398 ! delta_t_min = 0.001 390 399 391 400 ! 1. - Save initial values, initialize tendencies, initialize output fields … … 398 407 !! deltatw0(i, k) = deltatw(i, k) 399 408 !! deltaqw0(i, k) = deltaqw(i, k) 400 !! t env(i, k) = tenv0(i, k)401 !! q e(i, k) = qe0(i, k)409 !! tb(i, k) = tb0(i, k) 410 !! qb(i, k) = qb0(i, k) 402 411 !! dtls(i, k) = 0. 403 412 !! dqls(i, k) = 0. 404 413 !! d_deltat_gw(i, k) = 0. 405 !! d_t env(i, k) = 0.406 !! d_q e(i, k) = 0.414 !! d_tb(i, k) = 0. 415 !! d_qb(i, k) = 0. 407 416 !! d_deltatw(i, k) = 0. 408 417 !! d_deltaqw(i, k) = 0. … … 416 425 deltatw0(:,:) = deltatw(:,:) 417 426 deltaqw0(:,:) = deltaqw(:,:) 418 t env(:,:) = tenv0(:,:)419 q e(:,:) = qe0(:,:)427 tb(:,:) = tb0(:,:) 428 qb(:,:) = qb0(:,:) 420 429 dtls(:,:) = 0. 421 430 dqls(:,:) = 0. 422 431 d_deltat_gw(:,:) = 0. 423 d_t env(:,:) = 0.424 d_q e(:,:) = 0.432 d_tb(:,:) = 0. 433 d_qb(:,:) = 0. 425 434 d_deltatw(:,:) = 0. 426 435 d_deltaqw(:,:) = 0. … … 544 553 !<jyg 545 554 dth(:,:) = 0. 546 t u(:,:) = 0.547 q u(:,:) = 0.555 tx(:,:) = 0. 556 qx(:,:) = 0. 548 557 dtke(:,:) = 0. 549 558 dqke(:,:) = 0. … … 591 600 ktop(i) = 0 592 601 kupper(i) = 0 593 sum_th u(i) = 0.594 sum_t u(i) = 0.595 sum_q u(i) = 0.596 sum_thv u(i) = 0.602 sum_thx(i) = 0. 603 sum_tx(i) = 0. 604 sum_qx(i) = 0. 605 sum_thvx(i) = 0. 597 606 sum_dth(i) = 0. 598 607 sum_dq(i) = 0. 599 sum_rho(i) = 0.600 608 sum_dtdwn(i) = 0. 601 609 sum_dqdwn(i) = 0. 602 610 603 av_th u(i) = 0.604 av_t u(i) = 0.605 av_q u(i) = 0.606 av_thv u(i) = 0.611 av_thx(i) = 0. 612 av_tx(i) = 0. 613 av_qx(i) = 0. 614 av_thvx(i) = 0. 607 615 av_dth(i) = 0. 608 616 av_dq(i) = 0. 609 av_rho(i) = 0.610 617 av_dtdwn(i) = 0. 611 618 av_dqdwn(i) = 0. … … 620 627 DO k = 1, klev 621 628 DO i = 1, klon 622 ! write(*,*)'wake 1',i,k,RD,t env(i,k)623 rho(i, k) = p(i, k)/(RD*t env(i,k))629 ! write(*,*)'wake 1',i,k,RD,tb(i,k) 630 rho(i, k) = p(i, k)/(RD*tb(i,k)) 624 631 ! write(*,*)'wake 2',rho(i,k) 625 632 IF (k==1) THEN 626 ! write(*,*)'wake 3',i,k,rd,t env(i,k)627 rhoh(i, k) = ph(i, k)/(RD*t env(i,k))628 ! write(*,*)'wake 4',i,k,rd,t env(i,k)633 ! write(*,*)'wake 3',i,k,rd,tb(i,k) 634 rhoh(i, k) = ph(i, k)/(RD*tb(i,k)) 635 ! write(*,*)'wake 4',i,k,rd,tb(i,k) 629 636 zhh(i, k) = 0 630 637 ELSE 631 ! write(*,*)'wake 5',rd,(t env(i,k)+tenv(i,k-1))632 rhoh(i, k) = ph(i, k)*2./(RD*(t env(i,k)+tenv(i,k-1)))638 ! write(*,*)'wake 5',rd,(tb(i,k)+tb(i,k-1)) 639 rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1))) 633 640 ! write(*,*)'wake 6',(-rhoh(i,k)*RG)+zhh(i,k-1) 634 641 zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1) 635 642 END IF 636 643 ! write(*,*)'wake 7',ppi(i,k) 637 the(i, k) = tenv(i, k)/ppi(i, k) 638 thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) 639 tu(i, k) = tenv(i, k) - deltatw(i, k)*sigmaw(i) 640 qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) 641 ! write(*,*)'wake 8',(RD*(tenv(i,k)+deltatw(i,k))) 642 rhow(i, k) = p(i, k)/(RD*(tenv(i,k)+deltatw(i,k))) 644 thb(i, k) = tb(i, k)/ppi(i, k) 645 thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) 646 tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i) 647 qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i) 648 ! write(*,*)'wake 8',(RD*(tb(i,k)+deltatw(i,k))) 643 649 dth(i, k) = deltatw(i, k)/ppi(i, k) 644 650 END DO … … 650 656 n2(i, k) = 0 651 657 ELSE 652 n2(i, k) = amax1(0., -RG**2/th e(i,k)*rho(i,k)*(the(i,k+1)-the(i,k-1))/ &658 n2(i, k) = amax1(0., -RG**2/thb(i,k)*rho(i,k)*(thb(i,k+1)-thb(i,k-1))/ & 653 659 (p(i,k+1)-p(i,k-1))) 654 660 END IF … … 667 673 END DO 668 674 669 ! Calcul de la masse volumique moyenne de la colonne (bdlmd)670 ! -----------------------------------------------------------------671 672 DO k = 1, klev673 DO i = 1, klon674 epaisseur1(i, k) = 0.675 epaisseur2(i, k) = 0.676 END DO677 END DO678 679 DO i = 1, klon680 epaisseur1(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*RG) + 1.681 epaisseur2(i, 1) = -(ph(i,2)-ph(i,1))/(rho(i,1)*RG) + 1.682 rhow_moyen(i, 1) = rhow(i, 1)683 END DO684 685 DO k = 2, klev686 DO i = 1, klon687 epaisseur1(i, k) = -(ph(i,k+1)-ph(i,k))/(rho(i,k)*RG) + 1.688 epaisseur2(i, k) = epaisseur2(i, k-1) + epaisseur1(i, k)689 rhow_moyen(i, k) = (rhow_moyen(i,k-1)*epaisseur2(i,k-1)+rhow(i,k)* &690 epaisseur1(i,k))/epaisseur2(i, k)691 END DO692 END DO693 694 675 695 676 ! Choose an integration bound well above wake top … … 698 679 ! Determine Wake top pressure (Ptop) from buoyancy integral 699 680 ! -------------------------------------------------------- 681 700 682 Do i=1, klon 701 683 wk_adv(i) = .True. 702 684 Enddo 703 685 Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, & 704 dth, hw0, rho, &705 ktop, wk_adv )686 dth, hw0, rho, delta_t_min, & 687 ktop, wk_adv, h_zzz, ptop1, ktop1) 706 688 707 689 IF (prt_level>=10) THEN … … 736 718 ! -------------------------------------- 737 719 738 ! Initialize sum_thv uto 1st level virt. pot. temp.720 ! Initialize sum_thvx to 1st level virt. pot. temp. 739 721 740 722 DO i = 1, klon 741 723 z(i) = 1. 742 724 dz(i) = 1. 743 sum_thv u(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)725 sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i) 744 726 sum_dth(i) = 0. 745 727 END DO … … 751 733 ! LJYF : ecriture pas sympa avec un tableau z(i) qui n'est pas utilise come tableau 752 734 z(i) = z(i) + dz(i) 753 sum_th u(i) = sum_thu(i) + thu(i, k)*dz(i)754 sum_t u(i) = sum_tu(i) + tu(i, k)*dz(i)755 sum_q u(i) = sum_qu(i) + qu(i, k)*dz(i)756 sum_thv u(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)735 sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i) 736 sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i) 737 sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i) 738 sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i) 757 739 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 758 740 sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) 759 sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)760 741 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) 761 742 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) … … 777 758 778 759 DO i = 1, klon 779 av_th u(i) = sum_thu(i)/hw0(i)780 av_t u(i) = sum_tu(i)/hw0(i)781 av_q u(i) = sum_qu(i)/hw0(i)782 av_thv u(i) = sum_thvu(i)/hw0(i)760 av_thx(i) = sum_thx(i)/hw0(i) 761 av_tx(i) = sum_tx(i)/hw0(i) 762 av_qx(i) = sum_qx(i)/hw0(i) 763 av_thvx(i) = sum_thvx(i)/hw0(i) 783 764 ! av_thve = sum_thve/hw0 784 765 av_dth(i) = sum_dth(i)/hw0(i) 785 766 av_dq(i) = sum_dq(i)/hw0(i) 786 av_rho(i) = sum_rho(i)/hw0(i)787 767 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 788 768 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 789 769 790 770 wape(i) = -RG*hw0(i)*(av_dth(i)+ & 791 epsim1*(av_th u(i)*av_dq(i)+av_dth(i)*av_qu(i)+av_dth(i)*av_dq(i)))/av_thvu(i)771 epsim1*(av_thx(i)*av_dq(i)+av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i) 792 772 793 773 END DO … … 852 832 ! -------------------------- 853 833 DO i = 1, klon 854 q0_min(i) = min((q e(i,1)-sigmaw(i)*deltaqw(i,1)), &855 (q e(i,1)+(1.-sigmaw(i))*deltaqw(i,1)))834 q0_min(i) = min((qb(i,1)-sigmaw(i)*deltaqw(i,1)), & 835 (qb(i,1)+(1.-sigmaw(i))*deltaqw(i,1))) 856 836 END DO 857 837 DO k = 2, klev 858 838 DO i = 1, klon 859 q1_min(i) = min((q e(i,k)-sigmaw(i)*deltaqw(i,k)), &860 (q e(i,k)+(1.-sigmaw(i))*deltaqw(i,k)))839 q1_min(i) = min((qb(i,k)-sigmaw(i)*deltaqw(i,k)), & 840 (qb(i,k)+(1.-sigmaw(i))*deltaqw(i,k))) 861 841 IF (q1_min(i)<=q0_min(i)) THEN 862 842 q0_min(i) = q1_min(i) … … 890 870 ! ------------------------------------------------------------------------ 891 871 ! 892 wk_adv = .True.893 872 DO isubstep = 1, nsub 894 873 ! … … 896 875 ! wk_adv is the logical flag enabling wake evolution in the time advance 897 876 ! loop 898 899 877 DO i = 1, klon 900 878 wk_adv(i) = ok_qx_qw(i) .AND. alpha(i) >= 1. … … 1163 1141 IF (wk_adv(i)) THEN 1164 1142 omgtop(i) = -rho(i, ktop(i))*RG*omgtop(i) 1165 dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1)) 1143 !! LJYF dp_deltomg(i, 1) = omgtop(i)/(ptop(i)-ph(i,1)) 1144 dp_deltomg(i, 1) = omgtop(i)/min(ptop(i)-ph(i,1),-smallestreal) 1166 1145 END IF 1167 1146 END DO … … 1268 1247 IF (wk_adv(i) .AND. k<=kupper(i)+1) THEN 1269 1248 dth(i, k) = deltatw(i, k)/ppi(i, k) 1270 th1(i, k) = th e(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area1271 th2(i, k) = th e(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake1272 q1(i, k) = q e(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area1273 q2(i, k) = q e(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake1249 th1(i, k) = thb(i, k) - sigmaw(i)*dth(i, k) ! undisturbed area 1250 th2(i, k) = thb(i, k) + (1.-sigmaw(i))*dth(i, k) ! wake 1251 q1(i, k) = qb(i, k) - sigmaw(i)*deltaqw(i, k) ! undisturbed area 1252 q2(i, k) = qb(i, k) + (1.-sigmaw(i))*deltaqw(i, k) ! wake 1274 1253 END IF 1275 1254 END DO … … 1353 1332 ! C 1354 1333 ! ----------------------------------------------------------------- 1355 d_t env(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- &1334 d_tb(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)- & 1356 1335 rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_th2(i,k+1))/ & 1357 1336 (ph(i,k)-ph(i,k+1)) & … … 1359 1338 (ph(i,k)-ph(i,k+1)) )*ppi(i, k) 1360 1339 1361 d_q e(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- &1340 d_qb(i, k) = dtimesub*((rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)- & 1362 1341 rre2(i)*omg(i,k+1)*(1.-sigmaw(i))*d_q2(i,k+1))/ & 1363 1342 (ph(i,k)-ph(i,k+1)) & … … 1365 1344 (ph(i,k)-ph(i,k+1)) ) 1366 1345 ELSE IF (wk_adv(i) .AND. k==kupper(i)) THEN 1367 d_t env(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k)1368 1369 d_q e(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1)))1346 d_tb(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_th1(i,k)/(ph(i,k)-ph(i,k+1)))*ppi(i, k) 1347 1348 d_qb(i, k) = dtimesub*(rre1(i)*omg(i,k)*sigmaw(i)*d_q1(i,k)/(ph(i,k)-ph(i,k+1))) 1370 1349 1371 1350 END IF … … 1515 1494 ! Scale tendencies so that water vapour remains positive in w and x. 1516 1495 1517 CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, q e, d_qe, deltaqw, &1496 CALL wake_vec_modulation(klon, klev, wk_adv, epsilon_loc, qb, d_qb, deltaqw, & 1518 1497 d_deltaqw, sigmaw, d_sigmaw, alpha) 1519 1498 ! … … 1534 1513 DO i = 1, klon 1535 1514 IF (wk_adv(i) .AND. k<=kupper(i)) THEN 1536 d_t env(i, k) = alpha(i)*d_tenv(i, k)1537 d_q e(i, k) = alpha(i)*d_qe(i, k)1515 d_tb(i, k) = alpha(i)*d_tb(i, k) 1516 d_qb(i, k) = alpha(i)*d_qb(i, k) 1538 1517 d_deltatw(i, k) = alpha(i)*d_deltatw(i, k) 1539 1518 d_deltaqw(i, k) = alpha(i)*d_deltaqw(i, k) … … 1554 1533 DO i = 1, klon 1555 1534 IF (wk_adv(i) .AND. k<=kupper(i)) THEN 1556 dtls(i, k) = dtls(i, k) + d_t env(i, k)1557 dqls(i, k) = dqls(i, k) + d_q e(i, k)1535 dtls(i, k) = dtls(i, k) + d_tb(i, k) 1536 dqls(i, k) = dqls(i, k) + d_qb(i, k) 1558 1537 ! cc nrlmd 1559 1538 d_deltatw2(i, k) = d_deltatw2(i, k) + d_deltatw(i, k) … … 1566 1545 DO i = 1, klon 1567 1546 IF (wk_adv(i) .AND. k<=kupper(i)) THEN 1568 t env(i, k) = tenv0(i, k) + dtls(i, k)1569 q e(i, k) = qe0(i, k) + dqls(i, k)1570 th e(i, k) = tenv(i, k)/ppi(i, k)1547 tb(i, k) = tb0(i, k) + dtls(i, k) 1548 qb(i, k) = qb0(i, k) + dqls(i, k) 1549 thb(i, k) = tb(i, k)/ppi(i, k) 1571 1550 deltatw(i, k) = deltatw(i, k) + d_deltatw(i, k) 1572 1551 deltaqw(i, k) = deltaqw(i, k) + d_deltaqw(i, k) 1573 1552 dth(i, k) = deltatw(i, k)/ppi(i, k) 1574 ! c print*,'k,qx,qw',k,q e(i,k)-sigmaw(i)*deltaqw(i,k)1575 ! c $ ,q e(i,k)+(1-sigmaw(i))*deltaqw(i,k)1553 ! c print*,'k,qx,qw',k,qb(i,k)-sigmaw(i)*deltaqw(i,k) 1554 ! c $ ,qb(i,k)+(1-sigmaw(i))*deltaqw(i,k) 1576 1555 END IF 1577 1556 END DO … … 1732 1711 1733 1712 Call pkupper (klon, klev, ptop, ph, p, pupper, kupper, & 1734 dth, hw, rho, & 1735 ktop, wk_adv) 1736 1713 dth, hw, rho, delta_t_min, & 1714 ktop, wk_adv, h_zzz, ptop1, ktop1) 1737 1715 1738 1716 ! 5/ Set deltatw & deltaqw to 0 above kupper … … 1753 1731 DO i = 1, klon 1754 1732 IF (wk_adv(i)) THEN !!! nrlmd 1755 sum_th u(i) = 0.1756 sum_t u(i) = 0.1757 sum_q u(i) = 0.1758 sum_thv u(i) = 0.1733 sum_thx(i) = 0. 1734 sum_tx(i) = 0. 1735 sum_qx(i) = 0. 1736 sum_thvx(i) = 0. 1759 1737 sum_dth(i) = 0. 1760 1738 sum_dq(i) = 0. 1761 sum_rho(i) = 0.1762 1739 sum_dtdwn(i) = 0. 1763 1740 sum_dqdwn(i) = 0. 1764 1741 1765 av_th u(i) = 0.1766 av_t u(i) = 0.1767 av_q u(i) = 0.1768 av_thv u(i) = 0.1742 av_thx(i) = 0. 1743 av_tx(i) = 0. 1744 av_qx(i) = 0. 1745 av_thvx(i) = 0. 1769 1746 av_dth(i) = 0. 1770 1747 av_dq(i) = 0. 1771 av_rho(i) = 0.1772 1748 av_dtdwn(i) = 0. 1773 1749 av_dqdwn(i) = 0. … … 1778 1754 ! -------------------------------------- 1779 1755 1780 ! Initialize sum_thv uto 1st level virt. pot. temp.1756 ! Initialize sum_thvx to 1st level virt. pot. temp. 1781 1757 1782 1758 DO i = 1, klon … … 1784 1760 z(i) = 1. 1785 1761 dz(i) = 1. 1786 sum_thv u(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)1762 sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i) 1787 1763 sum_dth(i) = 0. 1788 1764 END IF … … 1795 1771 IF (dz(i)>0) THEN 1796 1772 z(i) = z(i) + dz(i) 1797 sum_th u(i) = sum_thu(i) + thu(i, k)*dz(i)1798 sum_t u(i) = sum_tu(i) + tu(i, k)*dz(i)1799 sum_q u(i) = sum_qu(i) + qu(i, k)*dz(i)1800 sum_thv u(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)1773 sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i) 1774 sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i) 1775 sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i) 1776 sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i) 1801 1777 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 1802 1778 sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) 1803 sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)1804 1779 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) 1805 1780 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) … … 1825 1800 DO i = 1, klon 1826 1801 IF (wk_adv(i)) THEN !!! nrlmd 1827 av_th u(i) = sum_thu(i)/hw0(i)1828 av_t u(i) = sum_tu(i)/hw0(i)1829 av_q u(i) = sum_qu(i)/hw0(i)1830 av_thv u(i) = sum_thvu(i)/hw0(i)1802 av_thx(i) = sum_thx(i)/hw0(i) 1803 av_tx(i) = sum_tx(i)/hw0(i) 1804 av_qx(i) = sum_qx(i)/hw0(i) 1805 av_thvx(i) = sum_thvx(i)/hw0(i) 1831 1806 av_dth(i) = sum_dth(i)/hw0(i) 1832 1807 av_dq(i) = sum_dq(i)/hw0(i) 1833 av_rho(i) = sum_rho(i)/hw0(i)1834 1808 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 1835 1809 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 1836 1810 1837 wape(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_th u(i)*av_dq(i) + &1838 av_dth(i)*av_q u(i)+av_dth(i)*av_dq(i)))/av_thvu(i)1811 wape(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + & 1812 av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i) 1839 1813 END IF 1840 1814 END DO … … 1886 1860 ! 1887 1861 END DO ! isubstep end sub-timestep loop 1888 1889 1862 ! 1890 1863 ! ------------------------------------------------------------------------ … … 1915 1888 ! cc 1916 1889 z(i) = 0. 1917 sum_th u(i) = 0.1918 sum_t u(i) = 0.1919 sum_q u(i) = 0.1920 sum_thv u(i) = 0.1890 sum_thx(i) = 0. 1891 sum_tx(i) = 0. 1892 sum_qx(i) = 0. 1893 sum_thvx(i) = 0. 1921 1894 sum_dth(i) = 0. 1922 1895 sum_half_dth(i) = 0. 1923 1896 sum_dq(i) = 0. 1924 sum_rho(i) = 0.1925 1897 sum_dtdwn(i) = 0. 1926 1898 sum_dqdwn(i) = 0. 1927 1899 1928 av_th u(i) = 0.1929 av_t u(i) = 0.1930 av_q u(i) = 0.1931 av_thv u(i) = 0.1900 av_thx(i) = 0. 1901 av_tx(i) = 0. 1902 av_qx(i) = 0. 1903 av_thvx(i) = 0. 1932 1904 av_dth(i) = 0. 1933 1905 av_dq(i) = 0. 1934 av_rho(i) = 0.1935 1906 av_dtdwn(i) = 0. 1936 1907 av_dqdwn(i) = 0. … … 1947 1918 IF (ok_qx_qw(i)) THEN 1948 1919 ! cc 1949 rho(i, k) = p(i, k)/(RD*t env(i,k))1920 rho(i, k) = p(i, k)/(RD*tb(i,k)) 1950 1921 IF (k==1) THEN 1951 rhoh(i, k) = ph(i, k)/(RD*t env(i,k))1922 rhoh(i, k) = ph(i, k)/(RD*tb(i,k)) 1952 1923 zhh(i, k) = 0 1953 1924 ELSE 1954 rhoh(i, k) = ph(i, k)*2./(RD*(t env(i,k)+tenv(i,k-1)))1925 rhoh(i, k) = ph(i, k)*2./(RD*(tb(i,k)+tb(i,k-1))) 1955 1926 zhh(i, k) = (ph(i,k)-ph(i,k-1))/(-rhoh(i,k)*RG) + zhh(i, k-1) 1956 1927 END IF 1957 the(i, k) = tenv(i, k)/ppi(i, k) 1958 thu(i, k) = (tenv(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) 1959 tu(i, k) = tenv(i, k) - deltatw(i, k)*sigmaw(i) 1960 qu(i, k) = qe(i, k) - deltaqw(i, k)*sigmaw(i) 1961 rhow(i, k) = p(i, k)/(RD*(tenv(i,k)+deltatw(i,k))) 1928 thb(i, k) = tb(i, k)/ppi(i, k) 1929 thx(i, k) = (tb(i,k)-deltatw(i,k)*sigmaw(i))/ppi(i, k) 1930 tx(i, k) = tb(i, k) - deltatw(i, k)*sigmaw(i) 1931 qx(i, k) = qb(i, k) - deltaqw(i, k)*sigmaw(i) 1962 1932 dth(i, k) = deltatw(i, k)/ppi(i, k) 1963 1933 END IF … … 1968 1938 ! ----------------------------------------------------------- 1969 1939 1970 ! Initialize sum_thv uto 1st level virt. pot. temp.1940 ! Initialize sum_thvx to 1st level virt. pot. temp. 1971 1941 1972 1942 DO i = 1, klon … … 1977 1947 dz(i) = 1. 1978 1948 dz_half(i) = 1. 1979 sum_thv u(i) = thu(i, 1)*(1.+epsim1*qu(i,1))*dz(i)1949 sum_thvx(i) = thx(i, 1)*(1.+epsim1*qx(i,1))*dz(i) 1980 1950 sum_dth(i) = 0. 1981 1951 END IF … … 1991 1961 IF (dz(i)>0) THEN 1992 1962 z(i) = z(i) + dz(i) 1993 sum_th u(i) = sum_thu(i) + thu(i, k)*dz(i)1994 sum_t u(i) = sum_tu(i) + tu(i, k)*dz(i)1995 sum_q u(i) = sum_qu(i) + qu(i, k)*dz(i)1996 sum_thv u(i) = sum_thvu(i) + thu(i, k)*(1.+epsim1*qu(i,k))*dz(i)1963 sum_thx(i) = sum_thx(i) + thx(i, k)*dz(i) 1964 sum_tx(i) = sum_tx(i) + tx(i, k)*dz(i) 1965 sum_qx(i) = sum_qx(i) + qx(i, k)*dz(i) 1966 sum_thvx(i) = sum_thvx(i) + thx(i, k)*(1.+epsim1*qx(i,k))*dz(i) 1997 1967 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) 1998 1968 sum_dq(i) = sum_dq(i) + deltaqw(i, k)*dz(i) 1999 sum_rho(i) = sum_rho(i) + rhow(i, k)*dz(i)2000 1969 sum_dtdwn(i) = sum_dtdwn(i) + dtdwn(i, k)*dz(i) 2001 1970 sum_dqdwn(i) = sum_dqdwn(i) + dqdwn(i, k)*dz(i) … … 2027 1996 IF (ok_qx_qw(i)) THEN 2028 1997 ! cc 2029 av_th u(i) = sum_thu(i)/hw0(i)2030 av_t u(i) = sum_tu(i)/hw0(i)2031 av_q u(i) = sum_qu(i)/hw0(i)2032 av_thv u(i) = sum_thvu(i)/hw0(i)1998 av_thx(i) = sum_thx(i)/hw0(i) 1999 av_tx(i) = sum_tx(i)/hw0(i) 2000 av_qx(i) = sum_qx(i)/hw0(i) 2001 av_thvx(i) = sum_thvx(i)/hw0(i) 2033 2002 av_dth(i) = sum_dth(i)/hw0(i) 2034 2003 av_dq(i) = sum_dq(i)/hw0(i) 2035 av_rho(i) = sum_rho(i)/hw0(i)2036 2004 av_dtdwn(i) = sum_dtdwn(i)/hw0(i) 2037 2005 av_dqdwn(i) = sum_dqdwn(i)/hw0(i) 2038 2006 2039 wape2(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_th u(i)*av_dq(i) + &2040 av_dth(i)*av_q u(i)+av_dth(i)*av_dq(i)))/av_thvu(i)2007 wape2(i) = -RG*hw0(i)*(av_dth(i)+epsim1*(av_thx(i)*av_dq(i) + & 2008 av_dth(i)*av_qx(i)+av_dth(i)*av_dq(i)))/av_thvx(i) 2041 2009 END IF 2042 2010 END DO … … 2374 2342 END SUBROUTINE wake 2375 2343 2376 SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, q e, d_qe, deltaqw, &2344 SUBROUTINE wake_vec_modulation(nlon, nl, wk_adv, epsilon_loc, qb, d_qb, deltaqw, & 2377 2345 d_deltaqw, sigmaw, d_sigmaw, alpha) 2378 2346 ! ------------------------------------------------------ … … 2384 2352 2385 2353 ! Input 2386 REAL q e(nlon, nl), d_qe(nlon, nl)2354 REAL qb(nlon, nl), d_qb(nlon, nl) 2387 2355 REAL deltaqw(nlon, nl), d_deltaqw(nlon, nl) 2388 2356 REAL sigmaw(nlon), d_sigmaw(nlon) … … 2410 2378 DO i = 1, nlon 2411 2379 IF (wk_adv(i)) THEN 2412 x = q e(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qe(i, k) + &2380 x = qb(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + d_qb(i, k) + & 2413 2381 (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - d_sigmaw(i) * & 2414 2382 (deltaqw(i,k)+d_deltaqw(i,k)) 2415 2383 a = -d_sigmaw(i)*d_deltaqw(i, k) 2416 b = d_q e(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - &2384 b = d_qb(i, k) + (zeta(i,k)-sigmaw(i))*d_deltaqw(i, k) - & 2417 2385 deltaqw(i, k)*d_sigmaw(i) 2418 c = q e(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon_loc2386 c = qb(i, k) + (zeta(i,k)-sigmaw(i))*deltaqw(i, k) + epsilon_loc 2419 2387 discrim = b*b - 4.*a*c 2420 2388 ! print*, 'x, a, b, c, discrim', x, a, b, c, discrim … … 2448 2416 2449 2417 SUBROUTINE pkupper (klon, klev, ptop, ph, p, pupper, kupper, & 2450 dth, hw_, rho, &2451 ktop, wk_adv )2418 dth, hw_, rho, delta_t_min, & 2419 ktop, wk_adv, h_zzz, ptop1, ktop1) 2452 2420 2453 2421 USE lmdz_wake_ini , ONLY : wk_pupper 2454 2422 USE lmdz_wake_ini , ONLY : RG 2455 2423 USE lmdz_wake_ini , ONLY : hwmin 2424 USE lmdz_wake_ini , ONLY : RD 2425 USE lmdz_wake_ini , ONLY : smallestreal 2426 2456 2427 IMPLICIT NONE 2457 2428 2458 INTEGER, INTENT(IN) :: klon,klev 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 2429 INTEGER, INTENT(IN) :: klon,klev 2430 REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: ph, p 2431 REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: rho 2432 LOGICAL, DIMENSION (klon) , INTENT(IN) :: wk_adv 2433 REAL, DIMENSION (klon,klev+1) , INTENT(IN) :: dth 2434 REAL, INTENT(IN) :: delta_t_min 2435 2436 REAL, DIMENSION (klon) , INTENT(OUT) :: hw_ 2437 REAL, DIMENSION (klon) , INTENT(OUT) :: ptop 2438 INTEGER, DIMENSION (klon) , INTENT(OUT) :: Ktop 2439 REAL, DIMENSION (klon) , INTENT(OUT) :: pupper 2440 INTEGER, DIMENSION (klon) , INTENT(OUT) :: kupper 2441 REAL, DIMENSION (klon) , INTENT(OUT) :: h_zzz !! 2442 REAL, DIMENSION (klon) , INTENT(OUT) :: Ptop1 !! 2443 INTEGER, DIMENSION (klon) , INTENT(OUT) :: ktop1 !! 2444 2445 2446 !===================================== 2447 ! local variables 2448 !===================================== 2449 2469 2450 INTEGER :: i,k 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 !INTEGER, SAVE :: compte=0 2479 2451 REAL, DIMENSION (klon) :: dthmin 2452 REAL, DIMENSION (klon) :: ptop_provis,ptop_new 2453 INTEGER, DIMENSION (klon) :: k_ptop_provis 2454 REAL, DIMENSION (klon) :: z, dz 2455 REAL, DIMENSION (klon) :: sum_dth 2456 REAL, DIMENSION (klon) :: omega !! 2457 REAL, DIMENSION (klon,klev+1) :: int_dth !! 2458 REAL, DIMENSION (klon,klev+1) :: dzz !! 2459 REAL, DIMENSION (klon,klev+1) :: zzz !! 2460 REAL, DIMENSION (klon) :: frac_int_dth !! 2461 REAL :: epsil !! 2462 REAL :: ddd!! 2463 2464 LOGICAL :: new_ptop 2465 2466 INTEGER, SAVE :: ipas=0 2480 2467 ! LJYF : a priori z, dz sum_dth sont aussi des variables internes 2481 2468 ! Les eliminer apres verification convergence numerique 2482 2469 2483 !compte=compte+1 2484 !print*,'compte=',compte 2485 2486 delta_t_min = 0.2 2470 ! delta_t_min = 0.2 2487 2471 2472 new_ptop=.FALSE. 2473 2488 2474 ! Determine Ptop from buoyancy integral 2489 2475 ! --------------------------------------- 2490 2476 2491 2477 ! - 1/ Pressure of the level where dth changes sign. 2492 !print*,'WAKE LJYF' 2493 2494 DO i = 1, klon 2478 ipas=ipas+1 2479 2480 DO i = 1, klon 2481 IF (wk_adv(i)) THEN 2495 2482 ptop_provis(i) = ph(i, 1) 2496 END DO 2497 2483 k_ptop_provis(i) = 1 2484 END IF 2485 END DO 2486 2487 2488 2498 2489 DO k = 2, klev 2499 2490 DO i = 1, klon 2500 ! if (compte==92) then2501 ! print*,'debug xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'2502 ! print*,'debug i =',i2503 ! print*,'debug k =',k2504 ! print*,'debug wk_adv(i) =',wk_adv(i)2505 ! print*,'debug ptop_provis(i) =',ptop_provis(i)2506 ! print*,'debug ph(i,1) =',ph(i,1)2507 ! print*,'debug dth(i,k) =',dth(i,k)2508 ! print*,'debug delta_t_min =',delta_t_min2509 ! print*,'debug p(i,k-1) =',p(i,k-1)2510 ! print*,'debug dth(i,k-1) =',dth(i,k-1)2511 ! print*,'debug p(i,k) =',p(i,k)2512 ! endif2513 2491 IF (wk_adv(i) .AND. ptop_provis(i)==ph(i,1) .AND. & 2514 2492 ! LJYF changer : dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 2515 dth(i,k)> -delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN2493 dth(i,k)>=-delta_t_min .AND. dth(i,k-1)<-delta_t_min) THEN 2516 2494 ptop_provis(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 2517 2495 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 2496 k_ptop_provis(i) = k 2518 2497 END IF 2519 2498 END DO … … 2534 2513 IF (wk_adv(i)) THEN 2535 2514 dz(i) = -(amax1(ph(i,k+1),ptop_provis(i))-ph(i,k))/(rho(i,k)*RG) 2536 IF (dz(i)>0 ) THEN2515 IF (dz(i)>0.) THEN 2537 2516 z(i) = z(i) + dz(i) 2538 2517 sum_dth(i) = sum_dth(i) + dth(i, k)*dz(i) … … 2577 2556 2578 2557 DO i = 1, klon 2558 IF (wk_adv(i)) THEN 2579 2559 ptop_new(i) = ptop(i) 2560 END IF 2580 2561 END DO 2581 2562 … … 2588 2569 ptop_new(i) = ((dth(i,k)+delta_t_min)*p(i,k-1) - & 2589 2570 (dth(i,k-1)+delta_t_min)*p(i,k))/(dth(i,k)-dth(i,k-1)) 2590 END IF 2591 END DO 2592 END DO 2593 2594 2595 DO i = 1, klon 2571 ! k_ptop_provis(i) = k 2572 END IF 2573 END DO 2574 END DO 2575 2576 2577 DO i = 1, klon 2578 IF (wk_adv(i)) THEN 2596 2579 ptop(i) = ptop_new(i) 2580 END IF 2597 2581 END DO 2598 2582 … … 2608 2592 ! PRINT *, 'wake-3, ktop(igout), kupper(igout) ', ktop(igout), kupper(igout) 2609 2593 ! ENDIF 2594 2595 2596 ! ----------------------------------------------------------------------- 2597 ! nouveau calcul de hw et ptop 2598 ! ----------------------------------------------------------------------- 2599 if (new_ptop) then 2600 2601 epsil = 0.05 ! 5 pour cent 2602 ! epsil = 0.20 2603 2604 DO i = 1, klon 2605 IF (wk_adv(i)) THEN 2606 int_dth(i,1) = 0. 2607 END IF 2608 END DO 2609 2610 2611 2612 DO K = 2, klev+1 2613 Do i = 1, klon 2614 IF (wk_adv(i)) THEN 2615 if (k<=k_ptop_provis(i)) then 2616 ddd=dth(i,k-1)*(ph(i,k-1) - min(ptop_provis(i),ph(i,k))) 2617 else 2618 ddd=0. 2619 endif 2620 int_dth(i,k) = int_dth(i,k-1) + ddd 2621 END IF 2622 END DO 2623 END DO 2624 ! print*, 'xxx, int_dth', (k,int_dth(1,k),k=1,klev) 2625 ! print*, 'xxx, k_ptop_provis', k_ptop_provis(1) 2626 2627 2628 DO i=1,klon 2629 IF (wk_adv(i)) THEN 2630 frac_int_dth(i)=(1.-epsil)*int_dth(i,k_ptop_provis(i)) 2631 ENDIF 2632 ENDDO 2633 DO k = 1,klev 2634 DO i =1, klon 2635 ! print*,ipas,'yyy ',k,int_dth(i,k),frac_int_dth(i) 2636 IF (wk_adv(i) .AND. int_dth(i,k)>=frac_int_dth(i)) THEN 2637 ktop1(i) = min(k, k_ptop_provis(i)) 2638 !print*,ipas,'yyy ktop1= ',ktop1 2639 END if 2640 END DO 2641 END DO 2642 !print*, 'LAMINE' 2643 2644 DO i = 1, klon 2645 IF (wk_adv(i)) THEN 2646 !print*, ipas,'xxx1, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ',ktop1 2647 ddd=int_dth(i,ktop1(i)+1)-int_dth(i,ktop1(i)) 2648 if (ddd==0.) then 2649 omega(i)=0. 2650 else 2651 omega(i) = (frac_int_dth(i) - int_dth(i,ktop1(i)))/ddd 2652 endif 2653 print*,'OMEGA ',omega(i) 2654 END IF 2655 END DO 2656 2657 print*, 'xxx' 2658 DO i = 1, klon 2659 IF (wk_adv(i)) THEN 2660 ! print*, 'xxx, int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) ', & 2661 ! int_dth(i,ktop1(i)), frac_int_dth(i), int_dth(i,ktop1(i)+1) 2662 ! print*, 'xxx, omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) ', & 2663 !e omega(i), ph(i,ktop1(i)), ph(i,ktop1(i)+1) 2664 ptop1(i) = min((1 - omega(i))*ph(i,ktop1(i)) + omega(i)*ph(i,ktop1(i)+1), ph(i,1)) 2665 END IF 2666 END DO 2667 2668 DO i=1, klon 2669 IF (wk_adv(i)) THEN 2670 zzz(i, 1) = 0 2671 END IF 2672 END DO 2673 DO k = 1, klev 2674 DO i = 1, klon 2675 IF (wk_adv(i)) THEN 2676 dzz(i,k) = (ph(i,k) - ph(i,k+1))/(rho(i,k)*RG) 2677 zzz(i,k+1) = zzz(i,k) + dzz(i,k) 2678 END IF 2679 END DO 2680 END DO 2681 2682 DO i =1, klon 2683 IF (wk_adv(i)) THEN 2684 h_zzz(i) = max((1- omega(i))*zzz(i,ktop1(i)) + omega(i)*zzz(i,ktop1(i)+1), hwmin) 2685 END IF 2686 END DO 2687 2688 endif 2689 !---------- FIN nouveau calcul hw et ptop ------------------------------------- 2690 2610 2691 2611 2692 kupper = 0
Note: See TracChangeset
for help on using the changeset viewer.
