Changeset 1127
- Timestamp:
- Dec 18, 2013, 9:50:25 AM (11 years ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/README
r1124 r1127 1961 1961 == 10/12/2013 == FGG 1962 1962 - Improved 15um cooling routines, with also better handling of errors. 1963 1964 == 18/12/2013 == EM 1965 - Bug fix in thermcall_main_mars, layer thickness "zdz" must be recomputed 1966 in every do ig=1,ngrid loop. -
trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90
r1033 r1127 326 326 !------------------------------------------------------------------ 327 327 ! 328 entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0. 329 lmin=1 328 entr_star(:,:)=0. ; detr_star(:,:)=0. 329 alim_star(:,:)=0. ; alim_star_tot(:)=0. 330 lmin(:)=1 330 331 331 332 !----------------------------------------------------------------------------- … … 419 420 420 421 ! is the thermal plume still active ? 421 422 do ig=1,ngrid 422 423 activecell(ig)=activecell(ig) & 423 424 & .and. zw2(ig,l)>1.e-10 & 424 425 & .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 425 426 enddo 426 427 427 428 !--------------------------------------------------------------------------- … … 437 438 !--------------------------------------------------------------------------- 438 439 439 440 441 442 443 444 445 446 447 448 449 450 440 do ig=1,ngrid 441 if(activecell(ig)) then 442 ztva_est(ig,l)=ztla(ig,l-1) 443 444 zdz=zlev(ig,l+1)-zlev(ig,l) 445 zbuoy(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) 446 447 ! Estimated vertical velocity squared 448 ! (discretized version of equation 12 in paragraph 40 of paper) 449 450 if (((a1*zbuoy(ig,l)/w_est(ig,l)-b1) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then 451 w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*a1*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*b1 & 451 452 & -2.*(1.-omega)*zdz*w_est(ig,l)*ae*(a1*zbuoy(ig,l)/w_est(ig,l)-b1)**be) 452 453 454 455 456 457 458 endif459 enddo453 else 454 w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*a1inv*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*b1inv) 455 endif 456 if (w_est(ig,l+1).lt.0.) then 457 w_est(ig,l+1)=zw2(ig,l) 458 endif 459 endif ! of if(activecell(ig)) 460 enddo ! of do ig=1,ngrid 460 461 461 462 !------------------------------------------------- … … 463 464 !------------------------------------------------- 464 465 465 do ig=1,ngrid466 if (activecell(ig)) then466 do ig=1,ngrid 467 if (activecell(ig)) then 467 468 468 469 zw2m=w_est(ig,l+1) 469 470 if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then 470 zdz=zlev(ig,l+1)-zlev(ig,l) 471 472 if((a1*(zbuoy(ig,l)/zw2m)-b1).gt.0.) then 471 473 472 474 ! ND entrainment rate, see equation 16 of paper (paragraph 43) 473 475 474 entr_star(ig,l)=f_star(ig,l)*zdz* &476 entr_star(ig,l)=f_star(ig,l)*zdz* & 475 477 & MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be) 478 476 479 else 477 entr_star(ig,l)=0.480 entr_star(ig,l)=0. 478 481 endif 479 482 … … 486 489 ! ND detrainment rate, see paragraph 44 of paper 487 490 488 detr_star(ig,l) = f_star(ig,l)*zdz* & 489 & ad 491 detr_star(ig,l) = f_star(ig,l)*zdz*ad 490 492 491 493 endif 492 494 else 493 detr_star(ig,l)=f_star(ig,l)*zdz* &495 detr_star(ig,l)=f_star(ig,l)*zdz* & 494 496 & MAX(ad,bd*zbuoy(ig,l)/zw2m) 495 497 … … 499 501 ! maximum between the source entrainment rate and the estimated entrainment rate. 500 502 501 if (l.lt.lalim(ig)) then502 alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))503 entr_star(ig,l)=0.504 endif503 if (l.lt.lalim(ig)) then 504 alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) 505 entr_star(ig,l)=0. 506 endif 505 507 506 508 ! Compute the non-dimensional upward mass flux at layer l+1 507 509 ! using equation 11 of appendix 4.2 in paper 508 510 509 f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &511 f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & 510 512 & -detr_star(ig,l) 511 513 512 endif513 enddo514 endif ! of if (activecell(ig)) 515 enddo ! of do ig=1,ngrid 514 516 515 517 ! ----------------------------------------------------------------------------------- … … 532 534 ! ----------------------------------------------------------------------------------- 533 535 ! ----------------------------------------------------------------------------------- 534 DO tic=0,5 536 DO tic=0,5 ! internal convergence loop 535 537 ! ----------------------------------------------------------------------------------- 536 538 ! ----------------------------------------------------------------------------------- … … 546 548 & (alim_star(ig,l)+entr_star(ig,l))*ztv(ig,l)) & 547 549 & /(f_star(ig,l+1)+detr_star(ig,l)) 548 549 endif 550 endif 550 551 enddo 551 552 … … 553 554 activetmp(:)=activetmp(:).and.(abs(ztla(:,l)-ztva(:,l)).gt.0.01) 554 555 555 ! Compute new buoyancu and vertical velocity 556 do ig=1,ngrid 557 if (activetmp(ig)) then 556 ! Compute new buoyancy and vertical velocity 557 do ig=1,ngrid 558 zdz=zlev(ig,l+1)-zlev(ig,l) 559 if (activetmp(ig)) then 558 560 ztva(ig,l) = ztla(ig,l) 559 561 zbuoy(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) 560 562 561 563 ! (discretized version of equation 12 in paragraph 40 of paper) 562 if (((a1*zbuoy(ig,l)/zw2(ig,l)-b1) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then 563 zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1*zbuoy(ig,l)- & 564 & 2.*zdz*zw2(ig,l)*b1-2.*(1.-omega)*zdz*zw2(ig,l)*ae*(a1*zbuoy(ig,l)/zw2(ig,l)-b1)**be) 564 if (((a1*zbuoy(ig,l)/zw2(ig,l)-b1) .gt. 0.) .and. & 565 (zw2(ig,l) .ne. 0.) ) then 566 zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1*zbuoy(ig,l)- & 567 2.*zdz*zw2(ig,l)*b1-2.*(1.-omega)*zdz*zw2(ig,l)* & 568 ae*(a1*zbuoy(ig,l)/zw2(ig,l)-b1)**be) 565 569 else 566 zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1inv*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*b1inv) 570 zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1inv*zbuoy(ig,l) & 571 -2.*zdz*zw2(ig,l)*b1inv) 567 572 endif 568 endif573 endif 569 574 enddo 570 575 … … 577 582 578 583 zw2m=zw2(ig,l+1) 584 zdz=zlev(ig,l+1)-zlev(ig,l) 579 585 if(zw2m .gt. 0) then 580 if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then 581 entr_star(ig,l)=f_star(ig,l)*zdz* & 582 & MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be) 586 if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then 587 entr_star(ig,l)=f_star(ig,l)*zdz* & 588 & MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be) 589 else 590 entr_star(ig,l)=0. 591 endif 592 593 if(zbuoy(ig,l) .gt. 0.) then 594 if(l .lt. lalim(ig)) then 595 596 detr_star(ig,l)=0. 597 598 else 599 detr_star(ig,l) = f_star(ig,l)*zdz*ad 600 601 endif 602 else 603 detr_star(ig,l)=f_star(ig,l)*zdz* & 604 & MAX(ad,bd*zbuoy(ig,l)/zw2m) 605 606 endif 583 607 else 584 entr_star(ig,l)=0. 585 endif 586 587 if(zbuoy(ig,l) .gt. 0.) then 588 if(l .lt. lalim(ig)) then 589 590 detr_star(ig,l)=0. 591 592 else 593 detr_star(ig,l) = f_star(ig,l)*zdz* & 594 & ad 595 596 endif 597 else 598 detr_star(ig,l)=f_star(ig,l)*zdz* & 599 & MAX(ad,bd*zbuoy(ig,l)/zw2m) 600 601 endif 602 else 603 entr_star(ig,l)=0. 604 detr_star(ig,l)=0. 605 endif 608 entr_star(ig,l)=0. 609 detr_star(ig,l)=0. 610 endif ! of if(zw2m .gt. 0) 606 611 607 612 ! If we are still in the source layer, we define the source layer entr. rate (alim_star) as the … … 616 621 ! using equation 11 of appendix 4.2 in paper 617 622 618 f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) &623 f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & 619 624 & -detr_star(ig,l) 620 625 621 endif622 enddo 626 endif ! of if (activetmp(ig)) 627 enddo ! of do ig=1,ngrid 623 628 ! ----------------------------------------------------------------------------------- 624 629 ! ----------------------------------------------------------------------------------- 625 ENDDO ! of internal convergence loop 630 ENDDO ! of internal convergence loop DO tic=0,5 626 631 ! ----------------------------------------------------------------------------------- 627 632 ! ----------------------------------------------------------------------------------- … … 632 637 633 638 do ig=1,ngrid 634 635 IF (thermverbose) THEN639 if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then 640 IF (thermverbose) THEN 636 641 print*,'thermcell_plume, particular case in velocity profile' 637 ENDIF638 639 642 ENDIF 643 zw2(ig,l+1)=0. 644 endif 640 645 641 646 if (zw2(ig,l+1).lt.0.) then 642 647 zw2(ig,l+1)=0. 643 648 endif 644 649 wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) 645 650 646 651 if (wa_moy(ig,l+1).gt.wmaxa(ig)) then … … 653 658 !========================================================================= 654 659 ! END OF THE LOOP ON VERTICAL LEVELS 655 enddo 660 enddo ! of do l=2,nlayer-1 656 661 !========================================================================= 657 662 !========================================================================= … … 708 713 do ig=1,ngrid 709 714 if ( zw2(ig,nlayer) > 1.e-10 ) then 710 print*,' WARNING !!!!! W2 non-zero in last layer'715 print*,'thermcell_main_mars: WARNING !!!!! W2 non-zero in last layer for ig=',ig 711 716 lmax(ig)=nlayer 712 717 endif … … 763 768 print*,'thermals have reached last layer of the model' 764 769 print*,'this is not good !' 770 zlmax=nlayer 765 771 endif 766 767 772 ! alim_star_clos is the source profile used for closure. It consists of the 768 773 ! modified source profile in the source layers, and the entrainment profile
Note: See TracChangeset
for help on using the changeset viewer.