- Timestamp:
- Aug 2, 2024, 2:12:03 PM (3 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_main.F90
r5135 r5158 197 197 !PRINT*,'thermcell_main debut' 198 198 ! WRITE(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)' 199 doig = 1, ngrid199 DO ig = 1, ngrid 200 200 f0(ig) = max(f0(ig), 1.e-2) 201 201 zmax0(ig) = max(zmax0(ig), 40.) … … 204 204 205 205 IF (prt_level>=20) THEN 206 doig = 1, ngrid206 DO ig = 1, ngrid 207 207 PRINT*, 'th_main ig f0', ig, f0(ig) 208 208 enddo … … 256 256 ! contenu thermcell_env : enddo 257 257 258 dol = 1, nlay259 doig = 1, ngrid258 DO l = 1, nlay 259 DO ig = 1, ngrid 260 260 zl(ig, l) = 0. 261 261 zu(ig, l) = puwind(ig, l) … … 296 296 !----------------------------------------------------------------------- 297 297 298 dol = 2, nlay298 DO l = 2, nlay 299 299 zlev(:, l) = 0.5 * (pphi(:, l) + pphi(:, l - 1)) / RG 300 300 enddo 301 301 zlev(:, 1) = 0. 302 302 zlev(:, nlay + 1) = (2. * pphi(:, nlay) - pphi(:, nlay - 1)) / RG 303 dol = 1, nlay303 DO l = 1, nlay 304 304 zlay(:, l) = pphi(:, l) / RG 305 305 enddo 306 dol = 1, nlay306 DO l = 1, nlay 307 307 deltaz(:, l) = zlev(:, l + 1) - zlev(:, l) 308 308 enddo … … 315 315 IF (prt_level>=10) WRITE(lunout, *) 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)' 316 316 rhobarz(:, 1) = rho(:, 1) 317 dol = 2, nlay317 DO l = 2, nlay 318 318 rhobarz(:, l) = 0.5 * (rho(:, l) + rho(:, l - 1)) 319 319 enddo 320 dol = 1, nlay320 DO l = 1, nlay 321 321 masse(:, l) = (pplev(:, l) - pplev(:, l + 1)) / RG 322 322 enddo … … 457 457 ! Le probleme vient du fait que linter et lmix sont souvent egaux a 1. 458 458 wmax_tmp = 0. 459 dol = 1, nlay459 DO l = 1, nlay 460 460 wmax_tmp(:) = max(wmax_tmp(:), zw2(:, l)) 461 461 enddo … … 562 562 ! Calcul de la fraction de l'ascendance 563 563 !------------------------------------------------------------------ 564 doig = 1, ngrid564 DO ig = 1, ngrid 565 565 fraca(ig, 1) = 0. 566 566 fraca(ig, nlay + 1) = 0. 567 567 enddo 568 dol = 2, nlay569 doig = 1, ngrid568 DO l = 2, nlay 569 DO ig = 1, ngrid 570 570 IF (zw2(ig, l)>1.e-10) THEN 571 571 fraca(ig, l) = fm(ig, l) / (rhobarz(ig, l) * zw2(ig, l)) … … 592 592 593 593 ! Temperature potentielle liquide effectivement mélangée par les thermiques 594 doll = 1, nlay595 doig = 1, ngrid594 DO ll = 1, nlay 595 DO ig = 1, ngrid 596 596 zthl(ig, ll) = ptemp(ig, ll) / zpspsk(ig, ll) 597 597 enddo … … 600 600 zthl, zdthladj, zta, lev_out) 601 601 602 doll = 1, nlay603 doig = 1, ngrid602 DO ll = 1, nlay 603 DO ig = 1, ngrid 604 604 z_o(ig, ll) = p_o(ig, ll) 605 605 enddo … … 610 610 #ifdef ISO 611 611 ! C Risi: on utilise directement la meme routine 612 doixt=1,ntiso613 doll=1,nlay612 DO ixt=1,ntiso 613 DO ll=1,nlay 614 614 DO ig=1,ngrid 615 615 xtpo_tmp(ig,ll)=xtpo(ixt,ig,ll) … … 619 619 CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, & 620 620 xtpo_tmp,xtpdoadj_tmp,xtzo_tmp,lev_out) 621 doll=1,nlay621 DO ll=1,nlay 622 622 DO ig=1,ngrid 623 623 xtpdoadj(ixt,ig,ll)=xtpdoadj_tmp(ig,ll) … … 674 674 675 675 ! PRINT*,'13 OK convect8' 676 dol = 1, nlay677 doig = 1, ngrid676 DO l = 1, nlay 677 DO ig = 1, ngrid 678 678 pdtadj(ig, l) = zdthladj(ig, l) * zpspsk(ig, l) 679 679 enddo … … 690 690 ! calcul du niveau de condensation 691 691 ! initialisation 692 doig = 1, ngrid692 DO ig = 1, ngrid 693 693 nivcon(ig) = 0 694 694 zcon(ig) = 0. 695 695 enddo 696 696 !nouveau calcul 697 doig = 1, ngrid697 DO ig = 1, ngrid 698 698 ! WARNING !!! verifier que c'est bien ztemp_env qu'on veut là 699 699 CHI = ztemp_env(ig, 1) / (1669.0 - 122.0 * z_o(ig, 1) / zqsat(ig, 1) - ztemp_env(ig, 1)) … … 701 701 enddo 702 702 !IM do k=1,nlay 703 dok = 1, nlay - 1704 doig = 1, ngrid703 DO k = 1, nlay - 1 704 DO ig = 1, ngrid 705 705 IF ((pcon(ig)<=pplay(ig, k)) & 706 706 .AND.(pcon(ig)>pplay(ig, k + 1))) THEN … … 711 711 !IM 712 712 ierr = 0 713 doig = 1, ngrid713 DO ig = 1, ngrid 714 714 IF (pcon(ig)<=pplay(ig, nlay)) THEN 715 715 zcon2(ig) = zlay(ig, nlay) - (pcon(ig) - pplay(ig, nlay)) / (RG * rho(ig, nlay)) / 100. … … 723 723 724 724 IF (prt_level>=1) PRINT*, '14b OK convect8' 725 dok = nlay, 1, -1726 doig = 1, ngrid725 DO k = nlay, 1, -1 726 DO ig = 1, ngrid 727 727 IF (zqla(ig, k)>1e-10) THEN 728 728 nivcon(ig) = k … … 734 734 !calcul des moments 735 735 !initialisation 736 dol = 1, nlay737 doig = 1, ngrid736 DO l = 1, nlay 737 DO ig = 1, ngrid 738 738 q2(ig, l) = 0. 739 739 wth2(ig, l) = 0. … … 746 746 IF (prt_level>=10)WRITE(lunout, *) & 747 747 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10' 748 dol = 1, nlay749 doig = 1, ngrid748 DO l = 1, nlay 749 DO ig = 1, ngrid 750 750 zf = fraca(ig, l) 751 751 zf2 = zf / (1. - zf) … … 765 765 enddo 766 766 !calcul des flux: q, thetal et thetav 767 dol = 1, nlay768 doig = 1, ngrid767 DO l = 1, nlay 768 DO ig = 1, ngrid 769 769 wq(ig, l) = fraca(ig, l) * zw2(ig, l) * (zqta(ig, l) * 1000. - p_o(ig, l) * 1000.) 770 770 wthl(ig, l) = fraca(ig, l) * zw2(ig, l) * (ztla(ig, l) - zthl(ig, l)) … … 779 779 ratqsdiff(:, :) = 0. 780 780 781 dol = 1, nlay782 doig = 1, ngrid781 DO l = 1, nlay 782 DO ig = 1, ngrid 783 783 IF (l<=lalim(ig)) THEN 784 784 var = var + alim_star(ig, l) * zqta(ig, l) * 1000. … … 789 789 IF (prt_level>=1) PRINT*, '14f OK convect8' 790 790 791 dol = 1, nlay792 doig = 1, ngrid791 DO l = 1, nlay 792 DO ig = 1, ngrid 793 793 IF (l<=lalim(ig)) THEN 794 794 zf = fraca(ig, l) … … 800 800 801 801 IF (prt_level>=1) PRINT*, '14g OK convect8' 802 dol = 1, nlay803 doig = 1, ngrid802 DO l = 1, nlay 803 DO ig = 1, ngrid 804 804 ratqsdiff(ig, l) = sqrt(vardiff) / (p_o(ig, l) * 1000.) 805 805 enddo … … 837 837 838 838 ! test sur la hauteur des thermiques ... 839 doi = 1, ngrid839 DO i = 1, ngrid 840 840 !IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) THEN 841 841 IF (prt_level>=10) THEN 842 842 PRINT*, 'WARNING ', comment, ' au point ', i, ' K= ', long(i) 843 843 PRINT*, ' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F* W2' 844 dok = 1, nlay844 DO k = 1, nlay 845 845 WRITE(6, '(i3,7f10.3)') k, pplay(i, k), ztv(i, k), 1000 * p_o(i, k), ztva(i, k), 1000 * zqla(i, k), f_star(i, k), zw2(i, k) 846 846 enddo … … 896 896 897 897 ! calcul du detrainement 898 dok = 1, nlay898 DO k = 1, nlay 899 899 detr0(:, k) = fm0(:, k) - fm0(:, k + 1) + entr0(:, k) 900 900 masse0(:, k) = (pplev(:, k) - pplev(:, k + 1)) / RG … … 907 907 detr(:, 1) = 0.5 * detr0(:, 1) 908 908 fm(:, 1) = 0. 909 dok = 1, nlay - 1909 DO k = 1, nlay - 1 910 910 masse(:, k + 1) = 0.5 * (masse0(:, k) + masse0(:, k + 1)) 911 911 entr(:, k + 1) = 0.5 * (entr0(:, k) + entr0(:, k + 1)) … … 917 917 q(:, :) = therm_tke_max(:, :) 918 918 !!! nrlmd le 16/09/2010 919 doig = 1, ngrid919 DO ig = 1, ngrid 920 920 qa(ig, 1) = q(ig, 1) 921 921 enddo … … 923 923 924 924 IF (1==1) THEN 925 dok = 2, nlay926 doig = 1, ngrid925 DO k = 2, nlay 926 DO ig = 1, ngrid 927 927 IF ((fm(ig, k + 1) + detr(ig, k)) * ptimestep> & 928 928 1.e-5 * masse(ig, k)) THEN … … 943 943 ! Calcul du flux subsident 944 944 945 dok = 2, nlay946 doig = 1, ngrid945 DO k = 2, nlay 946 DO ig = 1, ngrid 947 947 wqd(ig, k) = fm(ig, k) * q(ig, k) 948 948 IF (wqd(ig, k)<0.) THEN … … 951 951 enddo 952 952 enddo 953 doig = 1, ngrid953 DO ig = 1, ngrid 954 954 wqd(ig, 1) = 0. 955 955 wqd(ig, nlay + 1) = 0. … … 957 957 958 958 ! Calcul des tendances 959 dok = 1, nlay960 doig = 1, ngrid959 DO k = 1, nlay 960 DO ig = 1, ngrid 961 961 q(ig, k) = q(ig, k) + (detr(ig, k) * qa(ig, k) - entr(ig, k) * q(ig, k) & 962 962 - wqd(ig, k) + wqd(ig, k + 1)) &
Note: See TracChangeset
for help on using the changeset viewer.
