Changeset 628


Ignore:
Timestamp:
Apr 20, 2012, 12:11:47 PM (13 years ago)
Author:
acolaitis
Message:

Correction to advection problems in thermals + made the thermals model faster by limiting the vertical extension of loops with the height reached by thermals.

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90

    r625 r628  
    5757      REAL d_v_ajs(ngridmx,nlayermx)
    5858      REAL fm_therm(ngridmx,nlayermx+1), entr_therm(ngridmx,nlayermx)
    59       REAL detr_therm(ngridmx,nlayermx)
     59      REAL detr_therm(ngridmx,nlayermx),detrmod(ngridmx,nlayermx)
    6060      REAL zw2(ngridmx,nlayermx+1)
    6161      REAL fraca(ngridmx,nlayermx+1),zfraca(ngridmx,nlayermx+1)
    … …  
    6767      LOGICAL qtransport_thermals,dtke_thermals
    6868      INTEGER l,ig,iq,ii(1)
     69      CHARACTER (LEN=20) modname
    6970
    7071!--------------------------------------------------------
    … …  
    8788      REAL zzw2(ngridmx,nlayermx+1)
    8889      REAL zmax(ngridmx)
     90      INTEGER ndt,zlmax
    8991
    9092!--------------------------------------------------------
    … …  
    225227! Initialization of intermediary variables
    226228
    227          zfm_therm(:,:)=0.
    228          zentr_therm(:,:)=0.
    229          zdetr_therm(:,:)=0.
    230          zheatFlux(:,:)=0.
    231          zheatFlux_down(:,:)=0.
    232          zbuoyancyOut(:,:)=0.
    233          zbuoyancyEst(:,:)=0.
     229!         zfm_therm(:,:)=0. !init is done inside
     230!         zentr_therm(:,:)=0.
     231!         zdetr_therm(:,:)=0.
     232!         zheatFlux(:,:)=0. 
     233!         zheatFlux_down(:,:)=0.
     234!         zbuoyancyOut(:,:)=0.
     235!         zbuoyancyEst(:,:)=0.
    234236         zzw2(:,:)=0.
    235237         zmax(:)=0.
    236238         lmax(:)=0
    237          d_t_the(:,:)=0.
    238          d_u_the(:,:)=0.
    239          d_v_the(:,:)=0.
     239!         d_t_the(:,:)=0. !init is done inside
     240
     241!         d_u_the(:,:)=0. !transported outside
     242!         d_v_the(:,:)=0.
    240243         dq2_the(:,:)=0.
    241          if (nqmx .ne. 0) then
    242             d_q_the(:,:,:)=0.
     244
     245         if (nqmx .ne. 0 .and. ico2 .ne. 0) then
     246            d_q_the(:,:,ico2)=0.
    243247         endif
    244248
    … …  
    256260
    257261            d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact
    258             d_u_the(:,:)=d_u_the(:,:)*ptimestep*fact
    259             d_v_the(:,:)=d_v_the(:,:)*ptimestep*fact
    260 !            dq2_the(:,:)=dq2_the(:,:)*fact
     262!            d_u_the(:,:)=d_u_the(:,:)*ptimestep*fact
     263!            d_v_the(:,:)=d_v_the(:,:)*ptimestep*fact
     264            dq2_the(:,:)=dq2_the(:,:)*fact
    261265            if (ico2 .ne. 0) then
    262266               d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*ptimestep*fact
    263267            endif
    264268
    265              zmaxth(:)=zmaxth(:)+zmax(:)*fact
    266              lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
     269            zmaxth(:)=zmaxth(:)+zmax(:)*fact
     270            lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
    267271            fm_therm(:,:)=fm_therm(:,:)  &
    268272     &      +zfm_therm(:,:)*fact
    … …  
    288292
    289293           d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
    290            d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
    291            d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
     294!           d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
     295!           d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
    292296            if (ico2 .ne. 0) then
    293297               d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
    … …  
    297301
    298302            zt(:,:) = zt(:,:) + d_t_the(:,:)
    299             zu(:,:) = zu(:,:) + d_u_the(:,:)
    300             zv(:,:) = zv(:,:) + d_v_the(:,:)
     303!            zu(:,:) = zu(:,:) + d_u_the(:,:)
     304!            zv(:,:) = zv(:,:) + d_v_the(:,:)
    301305            if (ico2 .ne. 0) then
    302306             pq_therm(:,:,ico2) = &
    … …  
    310314
    311315      lmax(:)=nint(lmax_real(:))
     316      zlmax=MAXVAL(lmax(:))+2
     317      if (zlmax .ge. nlayermx) then
     318        print*,'thermals have reached last layer of the model'
     319        print*,'this is not good !'
     320      endif
     321
    312322
    313323! Now that we have computed total entrainment and detrainment, we can
    … …  
    321331      enddo
    322332
    323 !      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
    324 !     &      ,fm_therm,entr_therm,  &
    325 !     &     masse,zu,d_u_ajs)
    326 !
    327 !      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
    328 !     &       ,fm_therm,entr_therm,  &
    329 !     &     masse,zv,d_v_ajs)
     333      detrmod(:,:)=0.
     334      do l=1,zlmax
     335         do ig=1,ngridmx
     336            detrmod(ig,l)=fm_therm(ig,l)-fm_therm(ig,l+1) &
     337     &      +entr_therm(ig,l)
     338            if (detrmod(ig,l).lt.0.) then
     339               entr_therm(ig,l)=entr_therm(ig,l)-detrmod(ig,l)
     340               detrmod(ig,l)=0.
     341            endif
     342         enddo
     343      enddo
     344      ndt=10
     345      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
     346     &      ,fm_therm,entr_therm,detrmod,  &
     347     &     masse,zu,d_u_ajs,ndt,zlmax)
     348
     349      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
     350     &       ,fm_therm,entr_therm,detrmod,  &
     351     &     masse,zv,d_v_ajs,ndt,zlmax)
    330352
    331353      if (nqmx .ne. 0.) then
    … …  
    333355      if (iq .ne. ico2) then
    334356      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
    335      &     ,fm_therm,entr_therm,  &
    336      &    masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq))
     357     &     ,fm_therm,entr_therm,detrmod,  &
     358     &    masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),ndt,zlmax)
    337359      endif
    338360      ENDDO
    … …  
    340362
    341363      if (dtke_thermals) then
     364      detrmod(:,:)=0.
     365      ndt=1
     366      do l=1,zlmax
     367         do ig=1,ngridmx
     368            detrmod(ig,l)=fm_therm(ig,l)-fm_therm(ig,l+1) &
     369     &      +entr_therm(ig,l)
     370            if (detrmod(ig,l).lt.0.) then
     371               entr_therm(ig,l)=entr_therm(ig,l)-detrmod(ig,l)
     372               detrmod(ig,l)=0.
     373            endif
     374         enddo
     375      enddo
    342376      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
    343      &     ,fm_therm,entr_therm,  &
    344      &    masse,q2_therm,dq2_therm)
     377     &     ,fm_therm,entr_therm,detrmod,  &
     378     &    masse,q2_therm,dq2_therm,ndt,zlmax)
    345379      endif
    346380
    … …  
    363397! **********************************************************************
    364398
    365            pdu_th(:,:)=d_u_ajs(:,:)/ptimestep
    366            pdv_th(:,:)=d_v_ajs(:,:)/ptimestep
     399      do l=1,zlmax
     400        pdu_th(:,l)=d_u_ajs(:,l)
     401        pdv_th(:,l)=d_v_ajs(:,l)
     402      enddo
    367403
    368404           if(qtransport_thermals) then
    … …  
    370406               do iq=1,nqmx
    371407                if (iq .ne. ico2) then
    372                   pdq_th(:,:,iq)=d_q_ajs(:,:,iq)
     408                  do l=1,zlmax
     409                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq)
     410                  enddo
    373411                else
    374                   pdq_th(:,:,iq)=d_q_ajs(:,:,iq)/ptimestep
     412                  do l=1,zlmax
     413                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq)/ptimestep
     414                  enddo
    375415                endif
    376416               enddo
    … …  
    387427           ENDIF
    388428
    389            pdt_th(:,:)=d_t_ajs(:,:)/ptimestep
     429           do l=1,zlmax
     430              pdt_th(:,l)=d_t_ajs(:,l)/ptimestep
     431           enddo
    390432
    391433

  • trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90

    r621 r628  
    1       subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm0,entr0,  &
    2      &    masse0,q_therm,dq_therm)
     1      subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm,entr,detr,  &
     2     &    masse0,q_therm,dq_therm,ndt,zlmax)
    33      implicit none
    44
    5 ! #include "iniprint.h"
    65!=======================================================================
    76!
    … …  
    1110!   Version modifiee pour prendre les downdrafts a la place de la
    1211!   subsidence compensatoire
     12!
     13!   Version with sub-timestep for Martian thin layers
     14!
    1315!=======================================================================
    1416
    … …  
    2022      INTEGER, INTENT(IN) :: ngrid,nlayer
    2123      REAL, INTENT(IN) :: ptimestep
    22       REAL, INTENT(IN) :: fm0(ngridmx,nlayermx+1)
    23       REAL, INTENT(IN) ::entr0(ngridmx,nlayermx)
     24      REAL, INTENT(IN) :: fm(ngridmx,nlayermx+1)
     25      REAL, INTENT(IN) :: entr(ngridmx,nlayermx)
     26      REAL, INTENT(IN) :: detr(ngridmx,nlayermx)
    2427      REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx)
    2528      REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)
     29      INTEGER, INTENT(IN) :: ndt
     30      INTEGER, INTENT(IN) :: zlmax
    2631
    2732! ============================ OUTPUTS ===========================
    … …  
    3338      REAL q(ngridmx,nlayermx)
    3439      REAL qa(ngridmx,nlayermx)
    35       INTEGER ig,k
    36       REAL detr(ngridmx,nlayermx),entr(ngridmx,nlayermx)
    37       REAL wqd(ngridmx,nlayermx+1)
    38       REAL zzm
     40      INTEGER ig,k,i
     41      REAL invflux0(ngridmx,nlayermx)
     42      REAL ztimestep
    3943
    4044! =========== Init ==============================================
    … …  
    4246      qa(:,:)=q_therm(:,:)
    4347      q(:,:)=q_therm(:,:)
    44       entr(:,:)=entr0(:,:)
    45       detr(:,:)=0.
    4648
    47       do k=1,nlayermx
    48          do ig=1,ngridmx
    49             detr(ig,k)=fm0(ig,k)-fm0(ig,k+1)+entr(ig,k)
    50             if (detr(ig,k).lt.0.) then
    51                entr(ig,k)=entr(ig,k)-detr(ig,k)
    52                detr(ig,k)=0.
    53             endif
    54          enddo
    55       enddo
     49! ====== Computing q ============================================
    5650
    57 ! =========== Updraft ============================================
     51      dq_therm(:,:)=0.
     52      ztimestep=ptimestep/real(ndt)
     53      invflux0(:,:)=ztimestep/masse0(:,:)     
     54
     55      do i=1,ndt
    5856
    5957      do ig=1,ngrid
    … …  
    6159      enddo
    6260
    63       do k=2,nlayermx
     61      do k=2,zlmax
    6462         do ig=1,ngridmx
    65             if ((fm0(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
     63            if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
    6664     &         1.e-5*masse0(ig,k)) then
    67          qa(ig,k)=(fm0(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
    68      &     /(fm0(ig,k+1)+detr(ig,k))
     65         qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
     66     &     /(fm(ig,k+1)+detr(ig,k))
    6967            else
    7068               qa(ig,k)=q(ig,k)
    … …  
    7371      enddo
    7472
    75 ! =========== Subsidence ==========================================
    76 
    77       do k=2,nlayermx
    78          do ig=1,ngridmx
    79 
    80 ! Schema avec advection sur plus qu'une maille.
    81             zzm=masse0(ig,k)/ptimestep
    82             if (fm0(ig,k)>zzm) then
    83                wqd(ig,k)=(zzm*q(ig,k)+(fm0(ig,k)-zzm)*q(ig,k+1))
    84             else
    85                wqd(ig,k)=fm0(ig,k)*q(ig,k)
    86             endif
    87          enddo
    88       enddo
    89       do ig=1,ngrid
    90          wqd(ig,1)=0.
    91          wqd(ig,nlayermx+1)=0.
     73      do k=1,zlmax
     74           q(:,k)=q(:,k)+         &
     75     & (detr(:,k)*qa(:,k)-entr(:,k)*q(:,k) &
     76     &    -fm(:,k)*q(:,k)+fm(:,k+1)*q(:,k+1))  &
     77     &               *invflux0(:,k)
    9278      enddo
    9379
     80      enddo !of do i=1,ndt
    9481
    95 ! ====== dq ======================================================
     82! ====== Derivative ==============================================
    9683
    97       dq_therm(:,:)=0.
    9884
    99          do k=1,nlayermx-1
    100            q(:,k)=q(:,k)+         &
    101      & (detr(:,k)*qa(:,k)-entr(:,k)*q(:,k) &
    102      &       -wqd(:,k)+wqd(:,k+1))  &
    103      &               *ptimestep/masse0(:,k)
    104          enddo
    105 
    106          do k=1,nlayermx-1
     85         do k=1,zlmax
    10786          dq_therm(:,k)=(q(:,k)-q_therm(:,k))/ptimestep
    10887         enddo
    10988
     89! ==============
     90
    11091      return
    11192      end
     93

  • trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90

    r624 r628  
    9999      REAL f(ngridmx)
    100100
     101      REAL detrmod(ngridmx,nlayermx)
     102
    101103      REAL teta_th_int(ngridmx,nlayermx)
    102104      REAL teta_env_int(ngridmx,nlayermx)
    103105      REAL teta_down_int(ngridmx,nlayermx)
    104106
    105       CHARACTER (LEN=20) :: modname
    106107      CHARACTER (LEN=80) :: abort_message
     108      INTEGER ndt
    107109
    108110! ============= PLUME VARIABLES ============
    … …  
    123125      REAL denom(ngridmx)
    124126      REAL zlevinter(ngridmx)
     127      INTEGER zlmax
    125128
    126129! =========================================
    … …  
    186189      detr(:,:)=0.
    187190      fm(:,:)=0.
    188       zu(:,:)=pu(:,:)
    189       zv(:,:)=pv(:,:)
     191!      zu(:,:)=pu(:,:)
     192!      zv(:,:)=pv(:,:)
    190193      zhc(:,:)=pt(:,:)/zpopsk(:,:)
     194      ndt=1
    191195
    192196! **********************************************************************
    … …  
    888892! ===========================================================================
    889893
     894      zlmax=MAXVAL(lmax(:))+2
     895      if (zlmax .ge. nlayermx) then
     896        print*,'thermals have reached last layer of the model'
     897        print*,'this is not good !'
     898      endif
     899
    890900! Choix de la fonction d'alimentation utilisee pour la fermeture.
    891901
    … …  
    968978!-------------------------------------------------------------------------
    969979
    970       do l=1,nlayermx
     980      do l=1,zlmax
    971981         entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
    972982         detr(:,l)=f(:)*detr_star(:,l)
    973983      enddo
    974984
    975       do l=1,nlayermx
     985      do l=1,zlmax
    976986         do ig=1,ngridmx
    977987            if (l.lt.lmax(ig)) then
    … …  
    10071017!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    10081018
    1009       do l=1,nlayermx
     1019      do l=1,zlmax
    10101020
    10111021         do ig=1,ngridmx
    … …  
    11561166!-----------------------------------------------------------------------
    11571167
    1158       do l=1,nlayermx-1
     1168      do l=1,zlmax
    11591169         do ig=1,ngridmx
    11601170            eee0=entr(ig,l)
    … …  
    13471357! Calcul de la fraction de l'ascendance
    13481358!------------------------------------------------------------------
    1349       do ig=1,ngridmx
    1350          fraca(ig,1)=0.
    1351          fraca(ig,nlayermx+1)=0.
    1352       enddo
    1353       do l=2,nlayermx
     1359      fraca(:,:)=0.
     1360      do l=2,zlmax
    13541361         do ig=1,ngridmx
    13551362            if (zw2(ig,l).gt.1.e-10) then
    … …  
    15201527
    15211528      else
    1522 
    1523       call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
    1524      &      ,fm,entr,  &
    1525      &     masse,zu,pduadj)
    1526 
    1527       call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
    1528      &       ,fm,entr,  &
    1529      &     masse,zv,pdvadj)
     1529!      detrmod(:,:)=0.
     1530!      do k=1,zlmax
     1531!         do ig=1,ngridmx
     1532!            detrmod(ig,k)=fm(ig,k)-fm(ig,k+1) &
     1533!     &      +entr(ig,k)
     1534!            if (detrmod(ig,k).lt.0.) then
     1535!               entr(ig,k)=entr(ig,k)-detrmod(ig,k)
     1536!               detrmod(ig,k)=0.
     1537!            endif
     1538!         enddo
     1539!      enddo
     1540!
     1541!
     1542!      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
     1543!     &      ,fm,entr,detrmod,  &
     1544!     &     masse,zu,pduadj,ndt,zlmax)
     1545!
     1546!      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
     1547!     &       ,fm,entr,detrmod,  &
     1548!     &     masse,zv,pdvadj,ndt,zlmax)
    15301549
    15311550      endif
    … …  
    15381557! The rest is transported outside the sub-timestep loop
    15391558
     1559      ratiom(:,:)=1.
     1560
    15401561      if (ico2.ne.0) then
    1541 !      if (nqmx .ne. 0.) then
     1562      detrmod(:,:)=0.
     1563      do k=1,zlmax
     1564         do ig=1,ngridmx
     1565            detrmod(ig,k)=fm(ig,k)-fm(ig,k+1) &
     1566     &      +entr(ig,k)
     1567            if (detrmod(ig,k).lt.0.) then
     1568               entr(ig,k)=entr(ig,k)-detrmod(ig,k)
     1569               detrmod(ig,k)=0.
     1570            endif
     1571         enddo
     1572      enddo
     1573
    15421574      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
    1543      &     ,fm,entr,  &
    1544      &    masse,pq(:,:,ico2),pdqadj(:,:,ico2))
    1545 !      endif
     1575     &     ,fm,entr,detrmod,  &
     1576     &    masse,pq(:,:,ico2),pdqadj(:,:,ico2),ndt,zlmax)
    15461577
    15471578! Compute the ratio between theta and theta_m
    15481579     
    1549        do l=1,nlayermx
     1580       do l=1,zlmax
    15501581          do ig=1,ngridmx
    15511582             ratiom(ig,l)=1./(A*(pq(ig,l,ico2)+pdqadj(ig,l,ico2)*ptimestep)+B)
    15521583          enddo
    15531584       enddo
    1554        else
    1555           ratiom(:,:)=1.
     1585
    15561586       endif
    1557 
    15581587
    15591588!------------------------------------------------------------------
    … …  
    15611590!------------------------------------------------------------------
    15621591
    1563       do l=1,nlayermx
     1592      pdtadj(:,:)=0.
     1593      do l=1,zlmax
    15641594         do ig=1,ngridmx
    15651595         pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
    … …  
    16011631       teta_down_int(ig,nlayermx)=teta_down_int(ig,nlayermx-1)
    16021632      enddo
    1603       do l=1,nlayermx
     1633        heatFlux(:,:)=0.
     1634        buoyancyOut(:,:)=0.
     1635        buoyancyEst(:,:)=0.
     1636        heatFlux_down(:,:)=0.
     1637      do l=1,zlmax
    16041638       do ig=1,ngridmx
    16051639        heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))
Note: See TracChangeset for help on using the changeset viewer.