Changeset 161 for trunk


Ignore:
Timestamp:
Jun 16, 2011, 3:48:00 PM (13 years ago)
Author:
acolaitis
Message:

================================================
======== IMPLEMENTATION OF THERMALS ============
================================================

Author: A. Colaitis (2011-06-16)

The main goal of this revision is to start including the thermals into the model
for development purposes. Users should not use the thermals yet, as
several major configuration changes still need to be done.

This version includes :

  • updraft and downdraft parametrizations
  • velocity in the thermal, including drag
  • plume height analysis
  • closure equation
  • updraft transport of heat, tracers and momentum
  • downdraft transport of heat

This model should not be used without upcoming developments, namely :

  • downdraft transport of tracers and momentum
  • updraft & downdraft transport of q2 (tke)
  • revision of vdif_kc to compute q2 for non-stratified cases

Thermals could also include in a later revision :

  • momentum loss during transport (horizontal drag)

Compilation of the thermals has been successfully tested on ifort, gfortran and pgf90

================================================
================================================

M libf/phymars/callkeys.h
M libf/phymars/inifis.F

Added new control flags to call the thermals :

  • calltherm (false by default) <- to call thermals
  • outptherm (false by default) <- to output thermal-related diagnostics (for dev purposes)

================================================

M libf/phymars/vdifc.F
------> added a temporary output for thermal-related diagnostics

M libf/phymars/testphys1d.F
------> added treatment for a initialization from a profile of neutral gas (ar)

-> will be transformed in a decaying tracer for thermal diagnostics

M libf/phymars/physiq.F
------> added a section to call the thermals

-> changed the call to convadj
-> added thermal-related outputs for diagnostics

M libf/phymars/convadj.F
------> takes now into account the height of thermals to execute convective adjustment

=> note : convective adjustment needs to be activated when using thermals, in case of a

second instable layer above the thermals

================================================

A libf/phymars/calltherm_interface.F90
------> Interface between physiq.F and the thermals

A libf/phymars/calltherm_mars.F90
------> Routine running the sub-timestep of the thermals

A libf/phymars/thermcell_main_mars.F90
------> Main thermals routine specific to Martian physics

A libf/phymars/thermcell_dqupdown.F90
------> Thermals subroutine computing transport of quantities by updrafts and downdrafts

A libf/phymars/thermcell.F90
------> Module including parameters from the Earth to Mars importation. Will disappear in future dev

================================================
================================================

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
5 added
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/callkeys.h

    r38 r161  
    1212     &   ,callg2d,linear,rayleigh,tracer,active,doubleq,submicron       &
    1313     &   ,lifting,callddevil,scavenging,sedimentation,activice,water    &
    14      &   ,caps,photochem
     14     &   ,caps,photochem,calltherm,outptherm
    1515     
    1616      COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
    … …  
    2323     &   ,callstats,calleofdump                                         &
    2424     &   ,callnirco2,callnlte,callthermos,callconduct,                  &
    25      &    calleuv,callmolvis,callmoldiff,thermochem,thermoswater
     25     &    calleuv,callmolvis,callmoldiff,thermochem,thermoswater        &
     26     &   ,calltherm,outptherm
    2627
    2728

  • trunk/LMDZ.MARS/libf/phymars/convadj.F

    r38 r161  
    11      subroutine convadj(ngrid,nlay,nq,ptimestep,
    2      &                   pplay,pplev,ppopsk,
     2     &                   pplay,pplev,ppopsk,lmax_th,
    33     &                   pu,pv,ph,pq,
    44     &                   pdufi,pdvfi,pdhfi,pdqfi,
    … …  
    2121!     Modif. 2005 by F. Forget.
    2222!     Modif. 2010 by R. Wordsworth
    23 !     
     23!     Modif. 2011 by A. Colaitis
     24!
    2425!==================================================================
    2526
    … …  
    3839!     ---------
    3940
    40       INTEGER,intent(in) :: ngrid,nlay
     41      INTEGER,intent(in) :: ngrid,nlay,lmax_th(ngrid)
    4142      REAL,intent(in) :: ptimestep
    4243      REAL,intent(in) :: ph(ngrid,nlay)
    … …  
    163164      DO l=2,nlay
    164165        DO ig=1,ngrid
    165           IF(zhc(ig,l).LT.zhc(ig,l-1)) vtest(ig)=.true.
    166         ENDDO
     166          IF((zhc(ig,l).LT.zhc(ig,l-1))                                 $
     167     $  .AND. (l .GT. lmax_th(ig))) vtest(ig)=.true.
     168       ENDDO
    167169      ENDDO
    168170
    … …  
    198200!     Test loop upwards on l2
    199201
    200           DO
    201             l2 = l2 + 1
    202             IF (l2 .GT. nlay) EXIT
    203             IF (zhc(i, l2) .LT. zhc(i, l2-1)) THEN
    204  
     202      DO
     203       l2 = l2 + 1
     204       IF (l2 .GT. nlay) EXIT
     205       IF ((zhc(i, l2) .LT. zhc(i, l2-1)).AND.(l2 .GT. lmax_th(i))) THEN
     206
    205207!     l2 is the highest level of the unstable column
    206208 
    … …  
    228230!     do we have to extend the column downwards?
    229231 
    230                 down = .false.
    231                 IF (l1 .ne. 1) then    !--  and then
    232                   IF (zhmc .lt. zhc(i, l1-1)) then
    233                     down = .true.
    234                   END IF
    235                 END IF
     232            down = .false.
     233            IF (l1 .ne. 1) then    !--  and then
     234              IF ((zhmc .lt. zhc(i, l1-1)).and.(l1.gt.lmax_th(i))) then
     235                down = .true.
     236              END IF
     237            END IF
    236238 
    237239                ! this could be a problem...

  • trunk/LMDZ.MARS/libf/phymars/inifis.F

    r79 r161  
    204204         write(*,*) " calldifv = ",calldifv
    205205
     206         write(*,*) "call thermals ?"
     207         calltherm=.false. ! default value
     208         call getin("calltherm",calltherm)
     209         write(*,*) " calltherm = ",calltherm
     210
     211         write(*,*) "output thermal diagnostics ?"
     212         outptherm=.false. ! default value
     213         call getin("outptherm",outptherm)
     214         write(*,*) " outptherm = ",outptherm
     215
    206216         write(*,*) "call convective adjustment ?"
    207217         calladj=.true. ! default value
    … …  
    209219         write(*,*) " calladj = ",calladj
    210220         
     221         if (calltherm .and. (.not. calladj)) then
     222          print*,'Convadj has to be activated when using thermals'
     223          stop
     224         endif
    211225
    212226         write(*,*) "call CO2 condensation ?"

  • trunk/LMDZ.MARS/libf/phymars/physiq.F

    r83 r161  
    66     $            pw,
    77     $            pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
    8 
    98
    109      IMPLICIT NONE
    … …  
    301300      REAL time_phys
    302301
     302c Variables from thermal
     303
     304      REAL lmax_th_out(ngrid)
     305      REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
     306      REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
     307      INTEGER lmax_th(ngrid)
     308      REAL dtke_th(ngrid,nlayer+1)
     309      REAL dummycol(ngrid)
    303310c=======================================================================
    304311
    … …  
    515522c          Radiative transfer
    516523c          ------------------
     524
    517525           CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
    518526     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
    … …  
    599607c    4. Vertical diffusion (turbulent mixing):
    600608c    -----------------------------------------
    601 c
     609
    602610      IF (calldifv) THEN
    603 
    604611
    605612         DO ig=1,ngrid
    … …  
    624631     &        zdqdif,zdqsdif)
    625632
     633
    626634         DO l=1,nlayer
    627635            DO ig=1,ngrid
    … …  
    661669      ENDIF ! of IF (calldifv)
    662670
     671c-----------------------------------------------------------------------
     672c   TEST. Thermals :
     673c HIGHLY EXPERIMENTAL, BEWARE !!
     674c   -----------------------------
     675 
     676      if(calltherm) then
     677 
     678        call calltherm_interface(ngrid,nlayer,firstcall,
     679     $ long,lati,zzlev,zzlay,
     680     $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
     681     $ pplay,pplev,pphi,nq,zpopsk,
     682     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,dtke_th)
     683 
     684         DO l=1,nlayer
     685           DO ig=1,ngrid
     686              pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l)
     687              pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l)
     688              pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l)
     689              q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep
     690           ENDDO
     691        ENDDO
     692 
     693        DO ig=1,ngrid
     694          q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep
     695        ENDDO     
     696   
     697        if (tracer) then
     698        DO iq=1,nq
     699         DO l=1,nlayer
     700           DO ig=1,ngrid
     701             pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq)
     702           ENDDO
     703         ENDDO
     704        ENDDO
     705        endif
     706
     707        else   !of if calltherm
     708        lmax_th(:)=0
     709        end if
    663710
    664711c-----------------------------------------------------------------------
    … …  
    678725
    679726         CALL convadj(ngrid,nlayer,nq,ptimestep,
    680      $                pplay,pplev,zpopsk,
     727     $                pplay,pplev,zpopsk,lmax_th,
    681728     $                pu,pv,zh,pq,
    682729     $                pdu,pdv,zdh,pdq,
    683730     $                zduadj,zdvadj,zdhadj,
    684731     $                zdqadj)
     732
    685733
    686734         DO l=1,nlayer
    … …  
    11441192c        which can later be used to make the statistic files of the run:
    11451193c        "stats")          only possible in 3D runs !
    1146 
    11471194         
    11481195         IF (callstats) THEN
    … …  
    12611308c        call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2,
    12621309c    &                  emis)
     1310         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay)
     1311         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev)
    12631312         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2,
    12641313     &                  tsurf)
    … …  
    12771326     &                  fluxtop_sw_tot)
    12781327         call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt)
    1279 c        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
    1280 c        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
    1281 c        call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
     1328        call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu)
     1329        call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv)
     1330        call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw)
    12821331         call WRITEDIAGFI(ngrid,"rho","density","none",3,rho)
    12831332c        call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2)
    1284 c        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
     1333        call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh)
    12851334c        call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay)
    12861335c        call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2,
    … …  
    12901339c        call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate',
    12911340c    &                   'w.m-2',3,zdtlw)
     1341c        CALL WRITEDIAGFI(ngridmx,'tauTESap',
     1342c     &                         'tau abs 825 cm-1',
     1343c     &                         '',2,tauTES)
     1344
    12921345
    12931346c        ----------------------------------------------------------
    … …  
    14021455
    14031456c        ----------------------------------------------------------
     1457c        Outputs of thermals
     1458c        ----------------------------------------------------------
     1459
     1460
     1461!        call WRITEDIAGFI(ngrid,'dtke',
     1462!     &              'tendance tke thermiques','m**2/s**2',
     1463!     &                         3,dtke_th)
     1464!        call WRITEDIAGFI(ngrid,'d_u_ajs',
     1465!     &              'tendance u thermiques','m/s',
     1466!     &                         3,pdu_th*ptimestep)
     1467!        call WRITEDIAGFI(ngrid,'d_v_ajs',
     1468!     &              'tendance v thermiques','m/s',
     1469!     &                         3,pdv_th*ptimestep)
     1470!        if (tracer) then
     1471!        if (nq .eq. 2) then
     1472!        call WRITEDIAGFI(ngrid,'deltaq_th',
     1473!     &              'delta q thermiques','kg/kg',
     1474!     &                         3,ptimestep*pdq_th(:,:,2))
     1475!        endif
     1476!        endif
     1477
     1478
     1479c        ----------------------------------------------------------
    14041480c        Output in netcdf file "diagsoil.nc" for subterranean
    14051481c          variables (output every "ecritphy", as for writediagfi)
    … …  
    14351511c        CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1')
    14361512
     1513! THERMALS STUFF 1D
     1514         if(calltherm) then
     1515
     1516        lmax_th_out(:)=real(lmax_th(:))
     1517
     1518      if (ngridmx .eq. 1) then
     1519        call WRITEDIAGFI(ngridmx,'lmax_th',
     1520     &              'hauteur du thermique','K',
     1521     &                         0,lmax_th_out)
     1522
     1523       else
     1524        call WRITEDIAGFI(ngridmx,'lmax_th',
     1525     &              'hauteur du thermique','K',
     1526     &                         2,lmax_th_out)
     1527
     1528       endif
     1529
     1530         co2col(:)=0.
     1531         dummycol(:)=0.
     1532         if (tracer) then
     1533         do l=1,nlayermx
     1534           do ig=1,ngrid
     1535             co2col(ig)=co2col(ig)+
     1536     &                  zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g
     1537         if (nqmx .gt. 1) then
     1538             dummycol(ig)=dummycol(ig)+
     1539     &                  zq(ig,l,2)*(pplev(ig,l)-pplev(ig,l+1))/g
     1540         endif
     1541         enddo
     1542         enddo
     1543
     1544         end if
     1545         call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass'          &
     1546     &                                      ,'kg/m-2',0,co2col)
     1547         call WRITEDIAGFI(ngrid,'dummycol','integrated dummy mass'      &
     1548     &                                      ,'kg/m-2',0,dummycol)
     1549         endif
     1550         call WRITEDIAGFI(ngrid,'w','vertical velocity'                 &
     1551     &                              ,'m/s',1,pw)
     1552         call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2)
     1553         call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0,
     1554     &                  tsurf)
     1555
     1556           call WRITEDIAGFI(ngrid,"fluxsurf_sw",
     1557     &                "Solar radiative flux to surface","W.m-2",0,
     1558     &                fluxsurf_sw_tot)
     1559
     1560         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
     1561         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)
    14371562! or output in diagfi.nc (for testphys1d)
    14381563         call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps)

  • trunk/LMDZ.MARS/libf/phymars/testphys1d.F

    r86 r161  
    229229            close(91)
    230230          endif ! of if (txt.eq."co2")
     231          ! Allow for an initial profile of argon
     232          ! Can also be used to introduce a decaying tracer
     233          ! in the 1D (TBD) to study thermals
     234          if (txt.eq."ar") then
     235            !look for a "profile_ar" input file
     236            open(91,file='profile_ar',status='old',
     237     &       form='formatted',iostat=ierr)
     238            if (ierr.eq.0) then
     239              read(91,*) qsurf(iq)
     240              do ilayer=1,nlayermx
     241                read(91,*) q(ilayer,iq)
     242              enddo
     243            else
     244              write(*,*) "No profile_ar file!"
     245            endif
     246            close(91)
     247          endif ! of if (txt.eq."ar")
     248
    231249          ! WATER VAPOUR
    232250          if (txt.eq."h2o_vap") then

  • trunk/LMDZ.MARS/libf/phymars/vdifc.F

    r77 r161  
    638638      ENDIF
    639639
     640      if (calltherm .and. outptherm) then
     641      if (ngrid .eq. 1) then
     642        call WRITEDIAGFI(ngrid,'zh','zh inside vdifc',
     643     &                       'SI',1,ph(:,:)+pdhfi(:,:)*ptimestep)
     644        call WRITEDIAGFI(ngrid,'zkh','zkh',
     645     &                       'SI',1,zkh)
     646      endif
     647      endif
     648
    640649      RETURN
    641650      END
Note: See TracChangeset for help on using the changeset viewer.