Changeset 1127 for trunk


Ignore:
Timestamp:
Dec 18, 2013, 9:50:25 AM (11 years ago)
Author:
emillour
Message:

Mars GCM:

  • Bug fix in thermcall_main_mars, layer thickness "zdz" must be recomputed in every do ig=1,ngrid loop.

EM

Location:
trunk/LMDZ.MARS
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r1124 r1127  
    19611961== 10/12/2013 == FGG
    19621962- 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  
    326326!------------------------------------------------------------------
    327327!
    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
    330331
    331332!-----------------------------------------------------------------------------
    … …  
    419420
    420421! is the thermal plume still active ?
    421           do ig=1,ngrid
     422        do ig=1,ngrid
    422423             activecell(ig)=activecell(ig) &
    423424      &                 .and. zw2(ig,l)>1.e-10 &
    424425      &                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
    425           enddo
     426        enddo
    426427
    427428!---------------------------------------------------------------------------
    … …  
    437438!---------------------------------------------------------------------------
    438439
    439           do ig=1,ngrid
    440              if(activecell(ig)) then
    441                 ztva_est(ig,l)=ztla(ig,l-1)
    442 
    443                 zdz=zlev(ig,l+1)-zlev(ig,l)
    444                 zbuoy(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    445 
    446                 ! Estimated vertical velocity squared
    447                 ! (discretized version of equation 12 in paragraph 40 of paper)
    448 
    449                 if (((a1*zbuoy(ig,l)/w_est(ig,l)-b1) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then
    450                 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 &
     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 &
    451452     & -2.*(1.-omega)*zdz*w_est(ig,l)*ae*(a1*zbuoy(ig,l)/w_est(ig,l)-b1)**be)
    452                 else
    453                 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)
    454                 endif
    455                 if (w_est(ig,l+1).lt.0.) then
    456                 w_est(ig,l+1)=zw2(ig,l)
    457                 endif
    458              endif
    459           enddo
     453            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
    460461
    461462!-------------------------------------------------
    … …  
    463464!-------------------------------------------------
    464465
    465       do ig=1,ngrid
    466         if (activecell(ig)) then
     466        do ig=1,ngrid
     467         if (activecell(ig)) then
    467468
    468469          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
    471473
    472474         ! ND entrainment rate, see equation 16 of paper (paragraph 43)
    473475
    474           entr_star(ig,l)=f_star(ig,l)*zdz*  &
     476            entr_star(ig,l)=f_star(ig,l)*zdz*  &
    475477        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
     478         
    476479          else
    477           entr_star(ig,l)=0.
     480            entr_star(ig,l)=0.
    478481          endif
    479482
    … …  
    486489         ! ND detrainment rate, see paragraph 44 of paper
    487490
    488                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
    489             &  ad
     491                detr_star(ig,l) = f_star(ig,l)*zdz*ad
    490492
    491493             endif
    492494          else
    493           detr_star(ig,l)=f_star(ig,l)*zdz*                        &
     495            detr_star(ig,l)=f_star(ig,l)*zdz*                        &
    494496                &     MAX(ad,bd*zbuoy(ig,l)/zw2m)
    495497
    … …  
    499501! maximum between the source entrainment rate and the estimated entrainment rate.
    500502
    501        if (l.lt.lalim(ig)) then
    502           alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
    503           entr_star(ig,l)=0.
    504        endif
     503          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
    505507
    506508! Compute the non-dimensional upward mass flux at layer l+1
    507509! using equation 11 of appendix 4.2 in paper
    508510
    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)  &
    510512     &              -detr_star(ig,l)
    511513
    512       endif
    513       enddo
     514          endif ! of if (activecell(ig))
     515        enddo ! of do ig=1,ngrid
    514516
    515517! -----------------------------------------------------------------------------------
    … …  
    532534! -----------------------------------------------------------------------------------
    533535! -----------------------------------------------------------------------------------
    534       DO tic=0,5  ! internal convergence loop
     536      DO tic=0,5 ! internal convergence loop
    535537! -----------------------------------------------------------------------------------
    536538! -----------------------------------------------------------------------------------
    … …  
    546548     &            (alim_star(ig,l)+entr_star(ig,l))*ztv(ig,l))  &
    547549     &            /(f_star(ig,l+1)+detr_star(ig,l))
    548 
    549         endif
     550       endif
    550551      enddo
    551552
    … …  
    553554      activetmp(:)=activetmp(:).and.(abs(ztla(:,l)-ztva(:,l)).gt.0.01)
    554555
    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
    558560           ztva(ig,l) = ztla(ig,l)
    559561           zbuoy(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    560562
    561563          ! (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)
    565569           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)
    567572           endif
    568       endif
     573        endif
    569574      enddo
    570575
    … …  
    577582
    578583          zw2m=zw2(ig,l+1)
     584          zdz=zlev(ig,l+1)-zlev(ig,l)
    579585          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
    583607          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)
    606611
    607612! If we are still in the source layer, we define the source layer entr. rate  (alim_star) as the
    … …  
    616621! using equation 11 of appendix 4.2 in paper
    617622
    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)  &
    619624     &              -detr_star(ig,l)
    620625
    621       endif
    622       enddo
     626        endif ! of if (activetmp(ig))
     627      enddo ! of do ig=1,ngrid
    623628! -----------------------------------------------------------------------------------
    624629! -----------------------------------------------------------------------------------
    625       ENDDO   ! of internal convergence loop
     630      ENDDO   ! of internal convergence loop DO tic=0,5
    626631! -----------------------------------------------------------------------------------
    627632! -----------------------------------------------------------------------------------
    … …  
    632637
    633638      do ig=1,ngrid
    634             if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
    635       IF (thermverbose) THEN
     639        if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
     640          IF (thermverbose) THEN
    636641           print*,'thermcell_plume, particular case in velocity profile'
    637       ENDIF
    638                 zw2(ig,l+1)=0.
    639             endif
     642          ENDIF
     643          zw2(ig,l+1)=0.
     644        endif
    640645
    641646        if (zw2(ig,l+1).lt.0.) then
    642647           zw2(ig,l+1)=0.
    643648        endif
    644            wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
     649        wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
    645650
    646651        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
    … …  
    653658!=========================================================================
    654659! END OF THE LOOP ON VERTICAL LEVELS
    655       enddo
     660      enddo ! of do l=2,nlayer-1
    656661!=========================================================================
    657662!=========================================================================
    … …  
    708713      do ig=1,ngrid
    709714      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
    711716          lmax(ig)=nlayer
    712717      endif
    … …  
    763768        print*,'thermals have reached last layer of the model'
    764769        print*,'this is not good !'
     770        zlmax=nlayer
    765771      endif
    766 
    767772! alim_star_clos is the source profile used for closure. It consists of the
    768773! modified source profile in the source layers, and the entrainment profile
Note: See TracChangeset for help on using the changeset viewer.