- Timestamp:
- Jul 24, 2024, 2:54:37 PM (4 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_thermcell_plume_6A.f90
r5105 r5116 102 102 103 103 ! Initialisations des variables r?elles 104 if (1==1) then104 if (1==1) THEN 105 105 ztva(:, :) = ztv(:, :) 106 106 ztva_est(:, :) = ztva(:, :) … … 173 173 ! Le panache va prendre au debut les caracteristiques de l'air contenu 174 174 ! dans cette couche. 175 if (active(ig)) then175 if (active(ig)) THEN 176 176 ztla(ig, 1) = zthl(ig, 1) 177 177 zqta(ig, 1) = po(ig, 1) … … 217 217 do ig = 1, ngrid 218 218 ! PRINT*,'active',active(ig),ig,l 219 if(active(ig)) then219 IF(active(ig)) THEN 220 220 zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig)) 221 221 ztva_est(ig, l) = ztla(ig, l - 1) * zpspsk(ig, l) + RLvCp * zqla_est(ig, l) … … 247 247 248 248 !-------------------------------------------------- 249 if (iflag_thermals_ed<8) then249 if (iflag_thermals_ed<8) THEN 250 250 !-------------------------------------------------- 251 251 !AJ052014: J'ai remplac?? la boucle do par un do while … … 284 284 ztv2 = atv2 * zlt + btv2 285 285 286 if (ztv(ig, lt)>ztv1.and.ztv(ig, lt)<ztv2) then 287 286 if (ztv(ig, lt)>ztv1.and.ztv(ig, lt)<ztv2) THEN 288 287 !-------------------------------------------------- 289 288 !AJ052014: D??calage de zinv qui est entre le haut … … 293 292 zinv = zltdwn + zdz3 * factinv 294 293 295 if (zlmeldwn>=zinv) then294 if (zlmeldwn>=zinv) THEN 296 295 ztv_est(ig, l) = atv2 * zlmel + btv2 297 296 zbuoyjam(ig, l) = fact_shell * RG * (ztva_est(ig, l) - ztv_est(ig, l)) / ztv_est(ig, l) & 298 297 + (1. - fact_shell) * zbuoy(ig, l) 299 elseif (zlmelup>=zinv) then298 elseif (zlmelup>=zinv) THEN 300 299 ztv_est2 = atv2 * 0.5 * (zlmelup + zinv) + btv2 301 300 ztv_est1 = atv1 * 0.5 * (zinv + zlmeldwn) + btv1 … … 312 311 endif 313 312 314 else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then 315 316 if (zlmeldwn>zltdwn) then 313 else ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) THEN 314 if (zlmeldwn>zltdwn) THEN 317 315 zbuoyjam(ig, l) = fact_shell * RG * ((ztva_est(ig, l) - & 318 316 ztv(ig, lt)) / ztv(ig, lt)) + (1. - fact_shell) * zbuoy(ig, l) … … 330 328 ! & po(ig,lt))/po(ig,lt)+((zdz2-lmel)/zdz3)*(zqta(ig,l-1)- & 331 329 ! & po(ig,lt-1))/po(ig,lt-1)) 332 endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) then 333 334 else ! if (iflag_thermals_ed.lt.8) then 330 endif ! if (ztv(ig,lt).gt.ztv1.and.ztv(ig,lt).lt.ztv2) THEN 331 else ! if (iflag_thermals_ed.lt.8) THEN 335 332 lt = l + 1 336 333 zlt = zlev(ig, lt) … … 349 346 ztv(ig, lt)) / ztv(ig, lt) + (1. - coefzlmel) * (ztva_est(ig, l) - & 350 347 ztv(ig, lt - 1)) / ztv(ig, lt - 1)) + 0. * zbuoy(ig, l) 351 endif ! if (iflag_thermals_ed.lt.8) then 352 348 endif ! if (iflag_thermals_ed.lt.8) THEN 353 349 !------------------------------------------------ 354 350 !AJAM:nouveau calcul de w? … … 377 373 !AJ052014: J'ai remplac? w_est(ig,l) par zw2(ig,l) 378 374 !-------------------------------------------------- 379 if (iflag_thermals_ed==8) then375 if (iflag_thermals_ed==8) THEN 380 376 ! Ancienne version 381 377 ! w_est(ig,l+1)=Max(0.0001,(zdz/zdzbis)*(exp(-zw2fact)* & … … 408 404 endif 409 405 410 if (iflag_thermals_ed<6) then406 if (iflag_thermals_ed<6) THEN 411 407 zalpha = f0(ig) * f_star(ig, l) / sqrt(w_est(ig, l + 1)) / rhobarz(ig, l) 412 408 ! fact_epsilon=0.0005/(zalpha+0.025)**0.5 … … 431 427 !-------------------------------------------------- 432 428 433 ! if (w_est(ig,l+1).lt.0.) then429 ! if (w_est(ig,l+1).lt.0.) THEN 434 430 ! w_est(ig,l+1)=zw2(ig,l) 435 431 ! w_est(ig,l+1)=0.0001 … … 445 441 446 442 do ig = 1, ngrid 447 if (active(ig)) then 448 443 if (active(ig)) THEN 449 444 ! zw2m=max(0.5*(w_est(ig,l)+w_est(ig,l+1)),0.1) 450 445 zw2m = w_est(ig, l + 1) … … 476 471 ! entrbis=entr_star(ig,l) 477 472 478 if (iflag_thermals_ed<6) then473 if (iflag_thermals_ed<6) THEN 479 474 fact_epsilon = 0.0002 / (zalpha + 0.1) 480 475 endif … … 513 508 ! En dessous de lalim, on prend le max de alim_star et entr_star pour 514 509 ! alim_star et 0 sinon 515 if (l<lalim(ig)) then510 if (l<lalim(ig)) THEN 516 511 alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l)) 517 512 entr_star(ig, l) = 0. 518 513 endif 519 ! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) then514 ! if (l.lt.lalim(ig).and.alim_star(ig,l)>alim_star(ig,l-1)) THEN 520 515 ! alim_star(ig,l)=entrbis 521 516 ! endif … … 536 531 activetmp(:) = active(:) .and. f_star(:, l + 1)>1.e-10 537 532 do ig = 1, ngrid 538 if (activetmp(ig)) then533 if (activetmp(ig)) THEN 539 534 Zsat = .FALSE. 540 535 ztla(ig, l) = (f_star(ig, l) * ztla(ig, l - 1) + & … … 551 546 CALL thermcell_qsat(ngrid, activetmp, pplev(:, l), ztemp, zqta(:, l), zqsatth(:, l)) 552 547 do ig = 1, ngrid 553 if (activetmp(ig)) then548 if (activetmp(ig)) THEN 554 549 ! on ecrit de maniere conservative (sat ou non) 555 550 ! T = Tl +Lv/Cp ql … … 582 577 ! & (zw2(ig,l)-zdw2)+zdw2)+(zdzbis-zdz)/zdzbis* & 583 578 ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2)) 584 if (iflag_thermals_ed==8) then579 if (iflag_thermals_ed==8) THEN 585 580 zw2(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (zw2(ig, l) - zdw2) + zdw2) 586 581 else … … 591 586 ! & (exp(-zw2factbis)*(zw2(ig,l-1)-zdw2bis)+zdw2bis)) 592 587 593 if (iflag_thermals_ed<6) then588 if (iflag_thermals_ed<6) THEN 594 589 zalpha = f0(ig) * f_star(ig, l) / sqrt(zw2(ig, l + 1)) / rhobarz(ig, l) 595 590 ! fact_epsilon=0.0005/(zalpha+0.025)**0.5 … … 620 615 nbpb = 0 621 616 do ig = 1, ngrid 622 if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) then617 if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) THEN 623 618 ! stop 'On tombe sur le cas particulier de thermcell_dry' 624 619 ! PRINT*,'On tombe sur le cas particulier de thermcell_plume' … … 628 623 endif 629 624 630 if (zw2(ig, l + 1)<0.) then625 if (zw2(ig, l + 1)<0.) THEN 631 626 linter(ig) = (l * (zw2(ig, l + 1) - zw2(ig, l)) & 632 627 - zw2(ig, l)) / (zw2(ig, l + 1) - zw2(ig, l)) 633 628 zw2(ig, l + 1) = 0. 634 629 !+CR:04/05/12:correction calcul linter pour calcul de zmax continu 635 elseif (f_star(ig, l + 1)<0.) then630 elseif (f_star(ig, l + 1)<0.) THEN 636 631 linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l)) & 637 632 - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l)) … … 642 637 wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1)) 643 638 644 if (wa_moy(ig, l + 1)>wmaxa(ig)) then639 if (wa_moy(ig, l + 1)>wmaxa(ig)) THEN 645 640 ! lmix est le niveau de la couche ou w (wa_moy) est maximum 646 641 !on rajoute le calcul de lmix_bis 647 if (zqla(ig, l)<1.e-10) then642 if (zqla(ig, l)<1.e-10) THEN 648 643 lmix_bis(ig) = l + 1 649 644 endif … … 653 648 enddo 654 649 655 if (nbpb>0) then650 if (nbpb>0) THEN 656 651 PRINT*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume' 657 652 endif … … 772 767 773 768 ! Initialisations des variables reeles 774 if (1==1) then769 if (1==1) THEN 775 770 ztva(:, :) = ztv(:, :) 776 771 ztva_est(:, :) = ztva(:, :) … … 824 819 do l = 1, nlay - 1 825 820 do ig = 1, ngrid 826 if (ztv(ig, l)> ztv(ig, l + 1) .and. ztv(ig, 1)>=ztv(ig, l)) then821 if (ztv(ig, l)> ztv(ig, l + 1) .and. ztv(ig, 1)>=ztv(ig, l)) THEN 827 822 alim_star(ig, l) = MAX((ztv(ig, l) - ztv(ig, l + 1)), 0.) & 828 823 * sqrt(zlev(ig, l + 1)) … … 834 829 do l = 1, nlay 835 830 do ig = 1, ngrid 836 if (alim_star_tot(ig) > 1.e-10) then831 if (alim_star_tot(ig) > 1.e-10) THEN 837 832 alim_star(ig, l) = alim_star(ig, l) / alim_star_tot(ig) 838 833 endif … … 853 848 ! Le panache va prendre au debut les caracteristiques de l'air contenu 854 849 ! dans cette couche. 855 if (active(ig)) then850 if (active(ig)) THEN 856 851 ztla(ig, 1) = zthl(ig, 1) 857 852 zqta(ig, 1) = po(ig, 1) … … 898 893 do ig = 1, ngrid 899 894 ! PRINT*,'active',active(ig),ig,l 900 if(active(ig)) then895 IF(active(ig)) THEN 901 896 zqla_est(ig, l) = max(0., zqta(ig, l - 1) - zqsat(ig)) 902 897 ztva_est(ig, l) = ztla(ig, l - 1) * zpspsk(ig, l) + RLvCp * zqla_est(ig, l) … … 916 911 w_est(ig, l + 1) = Max(0.0001, exp(-zw2fact) * (w_est(ig, l) - zdw2) + zdw2) 917 912 918 if (w_est(ig, l + 1)<0.) then913 if (w_est(ig, l + 1)<0.) THEN 919 914 w_est(ig, l + 1) = zw2(ig, l) 920 915 endif … … 928 923 929 924 do ig = 1, ngrid 930 if (active(ig)) then 931 925 if (active(ig)) THEN 932 926 zw2m = max(0.5 * (w_est(ig, l) + w_est(ig, l + 1)), 0.1) 933 927 zw2m = w_est(ig, l + 1) … … 948 942 ! En dessous de lalim, on prend le max de alim_star et entr_star pour 949 943 ! alim_star et 0 sinon 950 if (l<lalim(ig)) then944 if (l<lalim(ig)) THEN 951 945 alim_star(ig, l) = max(alim_star(ig, l), entr_star(ig, l)) 952 946 entr_star(ig, l) = 0. … … 966 960 activetmp(:) = active(:) .and. f_star(:, l + 1)>1.e-10 967 961 do ig = 1, ngrid 968 if (activetmp(ig)) then962 if (activetmp(ig)) THEN 969 963 Zsat = .FALSE. 970 964 ztla(ig, l) = (f_star(ig, l) * ztla(ig, l - 1) + & … … 982 976 983 977 do ig = 1, ngrid 984 if (activetmp(ig)) then978 if (activetmp(ig)) THEN 985 979 ! on ecrit de maniere conservative (sat ou non) 986 980 ! T = Tl +Lv/Cp ql … … 1010 1004 nbpb = 0 1011 1005 do ig = 1, ngrid 1012 if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) then1006 if (zw2(ig, l + 1)>0. .and. zw2(ig, l + 1)<1.e-10) THEN 1013 1007 ! stop 'On tombe sur le cas particulier de thermcell_dry' 1014 1008 ! PRINT*,'On tombe sur le cas particulier de thermcell_plume' … … 1018 1012 endif 1019 1013 1020 if (zw2(ig, l + 1)<0.) then1014 if (zw2(ig, l + 1)<0.) THEN 1021 1015 linter(ig) = (l * (zw2(ig, l + 1) - zw2(ig, l)) & 1022 1016 - zw2(ig, l)) / (zw2(ig, l + 1) - zw2(ig, l)) 1023 1017 zw2(ig, l + 1) = 0. 1024 elseif (f_star(ig, l + 1)<0.) then1018 elseif (f_star(ig, l + 1)<0.) THEN 1025 1019 linter(ig) = (l * (f_star(ig, l + 1) - f_star(ig, l)) & 1026 1020 - f_star(ig, l)) / (f_star(ig, l + 1) - f_star(ig, l)) … … 1031 1025 wa_moy(ig, l + 1) = sqrt(zw2(ig, l + 1)) 1032 1026 1033 if (wa_moy(ig, l + 1)>wmaxa(ig)) then1027 if (wa_moy(ig, l + 1)>wmaxa(ig)) THEN 1034 1028 ! lmix est le niveau de la couche ou w (wa_moy) est maximum 1035 1029 !on rajoute le calcul de lmix_bis 1036 if (zqla(ig, l)<1.e-10) then1030 if (zqla(ig, l)<1.e-10) THEN 1037 1031 lmix_bis(ig) = l + 1 1038 1032 endif … … 1042 1036 enddo 1043 1037 1044 if (nbpb>0) then1038 if (nbpb>0) THEN 1045 1039 PRINT*, 'WARNING on tombe ', nbpb, ' x sur un pb pour l=', l, ' dans thermcell_plume' 1046 1040 endif
Note: See TracChangeset
for help on using the changeset viewer.
