Ignore:
Timestamp:
Jan 26, 2012, 3:09:06 PM (13 years ago)
Author:
acolaitis
Message:

Added molar-mass vertical gradient computations in thermals, coupled to the advection scheme for co2 tracer of thermals. As a result, there are now gentle thermals in the polar night (associated with a wstar of about 1.5 m.s-1).

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r496 r508  
    2020#include "dimphys.h"
    2121#include "comcstfi.h"
     22#include "tracer.h"
     23#include "callkeys.h"
    2224
    2325! ============== INPUTS ==============
    … …  
    110112      REAL wmaxa(ngridmx)
    111113      REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
    112       LOGICAL active(ngridmx),activetmp(ngridmx)
     114      LOGICAL activecell(ngridmx),activetmp(ngridmx)
    113115      REAL a1,b1,ae,be,ad,bd,fdfu
    114116      INTEGER tic
    … …  
    163165! =========================================
    164166
     167! ============== Theta_M Variables ========
     168
     169      INTEGER ico2
     170      SAVE ico2
     171      REAL m_co2, m_noco2, A , B
     172      SAVE A, B
     173      LOGICAL firstcall
     174      save firstcall
     175      data firstcall/.true./
     176      REAL zhc(ngridmx,nlayermx)
     177      REAL zdzfull(ngridmx,nlayermx)
     178      REAL ratiom(ngridmx,nlayermx)
     179
     180! =========================================
     181
    165182!-----------------------------------------------------------------------
    166183!   initialisation:
    … …  
    172189!      zu(:,:)=pu(:,:)
    173190!      zv(:,:)=pv(:,:)
    174       ztv(:,:)=pt(:,:)/zpopsk(:,:)
     191      zhc(:,:)=pt(:,:)/zpopsk(:,:)
     192
     193! **********************************************************************
     194! Taking into account vertical molar mass gradients
     195! **********************************************************************
     196
     197      if(firstcall) then
     198        ico2=0
     199        if (tracer) then
     200!     Prepare Special treatment if one of the tracers is CO2 gas
     201           do iq=1,nqmx
     202             if (noms(iq).eq."co2") then
     203                ico2=iq
     204                m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)
     205                m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)
     206!               Compute A and B coefficient use to compute
     207!               mean molecular mass Mair defined by
     208!               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
     209!               1/Mair = A*q(ico2) + B
     210                A =(1/m_co2 - 1/m_noco2)
     211                B=1/m_noco2
     212             end if
     213           enddo
     214        endif
     215
     216      firstcall=.false.
     217      endif !of if firstcall
     218
     219!.......................................................................
     220!  Special treatment for co2
     221!.......................................................................
     222
     223      if (ico2.ne.0) then
     224!     Special case if one of the tracers is CO2 gas
     225         DO l=1,nlayermx
     226           DO ig=1,ngridmx
     227            ztv(ig,l) = zhc(ig,l)*(A*pq(ig,l,ico2)+B)
     228           ENDDO
     229         ENDDO
     230       else
     231          ztv(:,:)=zhc(:,:)
     232       end if
     233
    175234
    176235!------------------------------------------------------------------------
    … …  
    332391! couches sont instables.
    333392!-------------------------------------------------------------------------
    334       active(:)=ztv(:,1)>ztv(:,2)
     393      activecell(:)=ztv(:,1)>ztv(:,2)
    335394          do ig=1,ngridmx
    336395            if (ztv(ig,1)>=(ztv(ig,2))) then
    … …  
    380439! Le panache va prendre au debut les caracteristiques de l'air contenu
    381440! dans cette couche.
    382           if (active(ig)) then
     441          if (activecell(ig)) then
    383442          ztla(ig,1)=ztv(ig,1)
    384443!cr: attention, prise en compte de f*(1)=1 => AC : what ? f*(1) =0. ! (d'ou f*(2)=a*(1)
    … …  
    407466! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
    408467          do ig=1,ngridmx
    409              active(ig)=active(ig) &
     468             activecell(ig)=activecell(ig) &
    410469      &                 .and. zw2(ig,l)>1.e-10 &
    411470      &                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
    … …  
    423482
    424483          do ig=1,ngridmx
    425              if(active(ig)) then
     484             if(activecell(ig)) then
    426485
    427486!                if(l .lt. lalim(ig)) then
    … …  
    453512
    454513      do ig=1,ngridmx
    455         if (active(ig)) then
     514        if (activecell(ig)) then
    456515
    457516          zw2m=w_est(ig,l+1)
    … …  
    459518          if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
    460519          entr_star(ig,l)=f_star(ig,l)*zdz*  &
    461         &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
     520!        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
     521         &  MAX(0.,log(1. + 0.027*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
    462522          else
    463523          entr_star(ig,l)=0.
    … …  
    480540! new param from continuity eq with a fit on dfdz
    481541                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
    482             &  ad
     542!            &  ad
     543             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
    483544
    484545!               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)   !svn baseline
    … …  
    492553          else
    493554          detr_star(ig,l)=f_star(ig,l)*zdz*                        &
    494                 &     bd*zbuoy(ig,l)/zw2m
     555!                &     bd*zbuoy(ig,l)/zw2m
     556             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
     557
    495558
    496559!               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)      !svn baseline
    … …  
    529592
    530593      DO tic=0,1  ! internal convergence loop
    531       activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
     594      activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10
    532595      do ig=1,ngridmx
    533596       if (activetmp(ig)) then
    … …  
    895958         enddo
    896959
    897 !  Les "optimisations" du flux sont desactivees : moins de bidouilles
     960!  Les "optimisations" du flux sont desactivecelles : moins de bidouilles
    898961!  je considere que celles ci ne sont pas justifiees ou trop delicates
    899962!  pour MARS, d'apres des observations LES.
    … …  
    13321395
    13331396!------------------------------------------------------------------
     1397!  calcul du transport vertical de traceurs
     1398!------------------------------------------------------------------
     1399
     1400! We only transport co2 tracer because it is coupled to the scheme through theta_m
     1401! The rest is transported outside the sub-timestep loop
     1402
     1403      if (ico2.ne.0) then
     1404!      if (nqmx .ne. 0.) then
     1405      do l=1,nlayermx
     1406         zdzfull(:,l)=zlev(:,l+1)-zlev(:,l)
     1407      enddo
     1408
     1409      modname='tracer'
     1410      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
     1411     &     ,fm,entr,detr,  &
     1412     &    masse,pq(:,:,ico2),pdqadj(:,:,ico2),modname,zdzfull)
     1413!      endif
     1414
     1415! Compute the ratio between theta and theta_m
     1416     
     1417       do l=1,nlayermx
     1418          do ig=1,ngridmx
     1419             ratiom(ig,l)=1./(A*(pq(ig,l,ico2)+pdqadj(ig,l,ico2)*ptimestep)+B)
     1420          enddo
     1421       enddo
     1422       else
     1423          ratiom(:,:)=1.
     1424       endif
     1425
     1426
     1427!------------------------------------------------------------------
    13341428!  incrementation dt
    13351429!------------------------------------------------------------------
    … …  
    13371431      do l=1,nlayermx
    13381432         do ig=1,ngridmx
    1339            pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)
    1340 !           pdtadj(ig,l)=zdthladj(ig,l)*zpopsk(ig,l)
    1341          enddo
    1342       enddo
    1343 
    1344 !------------------------------------------------------------------
    1345 !  calcul du transport vertical de traceurs
    1346 !------------------------------------------------------------------
    1347 
    1348 !      if (nqmx .ne. 0.) then
    1349 !      modname='tracer'
    1350 !      DO iq=1,nqmx
    1351 !          call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
    1352 !          &      masse,pq(:,:,iq),pdqadj(:,:,iq),ztvd,fm_down,ztv,modname,lmax)
    1353 !
    1354 !      ENDDO
    1355 !      endif
     1433         pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
     1434         enddo
     1435      enddo
    13561436
    13571437!------------------------------------------------------------------
    … …  
    13781458      do l=1,nlayermx-1
    13791459       do ig=1,ngridmx
    1380         teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l))
    1381         teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l))
    1382         teta_env_int(ig,l)=0.5*(ztv(ig,l+1)+ztv(ig,l))
     1460        teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l))*ratiom(ig,l)
     1461        teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l))*ratiom(ig,l)
     1462        teta_env_int(ig,l)=0.5*(ztv(ig,l+1)+ztv(ig,l))*ratiom(ig,l)
    13831463       enddo
    13841464      enddo
    … …  
    13911471       do ig=1,ngridmx
    13921472        heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))
    1393 !        buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    1394 !        buoyancyEst(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     1473        buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
     1474        buoyancyEst(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    13951475        heatFlux_down(ig,l)=fm_down(ig,l)*(teta_down_int(ig,l)-teta_env_int(ig,l))/rhobarz(ig,l)
    13961476       enddo
Note: See TracChangeset for help on using the changeset viewer.