Ignore:
Timestamp:
Nov 23, 2011, 10:14:21 AM (14 years ago)
Author:
aslmd
Message:

LMDZ.MARS : NEW NLTE MODEL FROM GRANADA AMIGOS

23/11/11 == FGG + MALV

New parameterization of the NLTE 15 micron cooling. The old parameterization is kept as an option, including or not variable atomic oxygen concentration. A new flag is introduced in callphys.def, nltemodel, to select which parameterization wants to be used (new one, old one with variable [O], or old one with fixed [O], see below). Includes many new subroutines and commons in phymars. Some existing routines are also modified:

-physiq.F. Call to the new subroutine NLTE_leedat in first call. Call to nltecool modified to allow for variable atomic oxygen. Depending on the value of nltemodel, the new subroutine NLTEdlvr09_TCOOL is called instead of nltecool.

-inifis.F. Reading of nltemodel is added.

-callkeys.h Declaration of nltemodel is added.

The following lines should be added to callphys.def (ideally after setting callnlte):

# NLTE 15um scheme to use.
# 0-> Old scheme, static oxygen
# 1-> Old scheme, dynamic oxygen
# 2-> New scheme
nltemodel = 2

A new directory, NLTEDAT, has to be included in datagcm.

Improvements into NLTE NIR heating parameterization to take into account variability with O/CO2 ratio and SZA. A new subroutine, NIR_leedat.F, and a new common, NIRdata.h, are included. A new flag, nircorr, is added in callphys.def, to include or not these corrections. The following files are modified:

-nirco2abs.F: nq and pq are added in the arguments. The corrections factors are interpolated to the GCM grid and included in the clculation. A new subroutine, interpnir, is added at the end of the file.

-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs

-inifis: Reading new flag nircorr.

-callkeys.h: Declaration of nircorr.

The following lines have to be added to callphys.def (ideally after callnirco2):

# NIR NLTE correction for variable SZA and O/CO2?
# matters only if callnirco2=T
# 0-> no correction
# 1-> include correction
nircorr=1

A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.

Small changes to the molecular diffusion scheme to fix the number of species considered, to avoid problems when compiling with more than 15 tracers (for example, when CH4 is included). Modified subroutines: aeronomars/moldiff.F and aeronomars/moldiffcoeff.F

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

Legend:

Unmodified
Added
Removed

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

    r358 r414  
    1616     
    1717      COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
    18      &   ,dustbin,nqchem_min
     18     &   ,dustbin,nqchem_min,nltemodel,nircorr
    1919     
    2020      COMMON/callkeys_r/topdustref,solarcondate,semi,alphan
    … …  
    5252      logical photochem
    5353      integer nqchem_min
     54      integer nltemodel
     55      integer nircorr
    5456
    5557      integer swrtype ! type of short wave (solar wavelength) radiative

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

    r358 r414  
    134134
    135135         write(*,*) "Directory where external input files are:"
    136          datafile="/u/forget/WWW/datagcm/datafile"
     136         datafile="/data/fgglmd/datagcm/datafile"
    137137         call getin("datadir",datafile) ! default path
    138138         write(*,*) " datafile = ",trim(datafile)
    … …  
    223223         write(*,*) " callnlte = ",callnlte
    224224         
     225         nltemodel=0    !default value
     226         write(*,*) "NLTE model?"
     227         write(*,*) "0 -> old model, static O"
     228         write(*,*) "1 -> old model, dynamic O"
     229         write(*,*) "2 -> new model"
     230         write(*,*) "(matters only if callnlte=T)"
     231         call getin("nltemodel",nltemodel)
     232         write(*,*) " nltemodel = ",nltemodel
     233
    225234         write(*,*) "call CO2 NIR absorption ?",
    226235     &              "(matters only if callrad=T)"
    … …  
    228237         call getin("callnirco2",callnirco2)
    229238         write(*,*) " callnirco2 = ",callnirco2
    230          
     239
     240         write(*,*) "New NIR NLTE correction ?",
     241     $              "0-> old model (no correction)",
     242     $              "1-> new correction",
     243     $              "(matters only if callnirco2=T)"
     244         call getin("nircorr",nircorr)
     245         write(*,*) " nircorr = ",nircorr
     246
    231247         write(*,*) "call turbulent vertical diffusion ?"
    232248         calldifv=.true. ! default value

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

    r38 r414  
    1       SUBROUTINE nirco2abs(ngrid,nlayer,pplay,dist_sol,
     1      SUBROUTINE nirco2abs(ngrid,nlayer,pplay,dist_sol,nq,pq,
    22     $     mu0,fract,declin,pdtnirco2)
    33                                                   
    … …  
    1616c   see NLTE correction-factor of Lopez-Valverde et al (1998)
    1717c   Stephen Lewis 2000
    18 c 
     18c
     19c   jul 2011 malv+fgg    New corrections for NLTE implemented
    1920c   08/2002 : correction for bug when running with diurnal=F
    2021c
    … …  
    4849#include "callkeys.h"
    4950#include "comdiurn.h"
    50 
     51#include "NIRdata.h"
    5152
    5253c-----------------------------------------------------------------------
    … …  
    6263c    Local variables :
    6364c    -----------------
    64       INTEGER l,ig, n, nstep
     65      INTEGER l,ig, n, nstep,i
    6566      REAL co2heat0, zmu(ngridmx)
    6667
    … …  
    8485      parameter (n_p0=0.0015888279)
    8586
     87c     Variables added to implement NLTE correction factor (feb 2011)
     88      real    pyy(nlayer)
     89      real    cor1(nlayer),oldoco2(nlayer),alfa2(nlayer)
     90      real    p2011,cociente1,merge
     91      real    cor0,oco2gcm
     92      integer nq
     93      real    pq(ngrid,nlayer,nq)
     94
    8695c----------------------------------------------------------------------
    8796
    … …  
    97106         do ig=1,ngrid
    98107            zmu(ig)=sqrt(1224.*mu0(ig)*mu0(ig)+1.)/35.
    99          enddo
    100          do l=1,nlayer
    101            do ig=1,ngrid
    102              if(fract(ig).gt.0.) pdtnirco2(ig,l)=
     108
     109            if(nircorr.eq.1) then
     110               do l=1,nlayer
     111                  pyy(l)=pplay(ig,l)
     112               enddo
     113
     114               call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres)
     115               call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres)
     116               call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres)
     117            endif
     118
     119            do l=1,nlayer
     120!           Calculations for the O/CO2 correction
     121               if(nircorr.eq.1) then
     122                  cor0=1./(1.+n_p0/pplay(ig,l))**n_b
     123                  if(pq(ig,l,1).gt.1.e-6) then
     124                     oco2gcm=pq(ig,l,3)/pq(ig,l,1)
     125                  else
     126                     oco2gcm=1.e6
     127                  endif
     128                  cociente1=oco2gcm/oldoco2(l)
     129                  merge=alog10(cociente1)*alfa2(l)+alog10(cor0)*
     130     $                 (1.-alfa2(l))
     131                  merge=10**merge
     132                  p2011=sqrt(merge)*cor0
     133               else if (nircorr.eq.0) then
     134                  p2011=1.
     135                  cor1(l)=1.
     136               endif
     137
     138               if(fract(ig).gt.0.) pdtnirco2(ig,l)=
    103139     &             co2heat0*sqrt((700.*zmu(ig))/pplay(ig,l))
    104140     &             /(1.+n_p0/pplay(ig,l))**n_b
    105 
     141!           Corrections from tabulation
     142     $              * cor1(l) * p2011
    106143c           OLD SCHEME (forget et al. 1999)
    107144c    s           co2heat0*sqrt((700.*zmu(ig))/pplay(ig,l))
    … …  
    109146           enddo
    110147         enddo
     148         
    111149
    112150c     Averaging over diurnal cycle (if diurnal=F)
    … …  
    128166            do ig=1,ngrid
    129167               zmu(ig)=sqrt(1224.*mu0_int(ig)*mu0_int(ig)+1.)/35.
    130             enddo
    131             do l=1,nlayer
    132                do ig=1,ngrid
     168
     169               if(nircorr.eq.1) then
     170                  do l=1,nlayer
     171                     pyy(l)=pplay(ig,l)
     172                  enddo
     173                  call interpnir(cor1,pyy,nlayer,corgcm,pres1d,npres)
     174                  call interpnir(oldoco2,pyy,nlayer,oco21d,pres1d,npres)
     175                  call interpnir(alfa2,pyy,nlayer,alfa,pres1d,npres)
     176               endif
     177
     178               do l=1,nlayer
     179                  if(nircorr.eq.1) then
     180                     cor0=1./(1.+n_p0/pplay(ig,l))**n_b
     181                     oco2gcm=pq(ig,l,3)/pq(ig,l,1)
     182                     cociente1=oco2gcm/oldoco2(l)
     183                     merge=alog10(cociente1)*alfa2(l)+alog10(cor0)*
     184     $                    (1.-alfa2(l))
     185                     merge=10**merge
     186                     p2011=sqrt(merge)*cor0
     187                  else if (nircorr.eq.0) then
     188                     p2011=1.
     189                     cor1(l)=1.
     190                  endif
     191
    133192                  if(fract_int(ig).gt.0.) pdtnirco2(ig,l)=
    134193     &                 pdtnirco2(ig,l) + (1/float(nstep))*
    135194     &                 co2heat0*sqrt((700.*zmu(ig))/pplay(ig,l))
    136195     &                 /(1.+n_p0/pplay(ig,l))**n_b
     196!                      Corrections from tabulation
     197     $                 * cor1(l) * p2011
    137198               enddo
    138199            enddo
    … …  
    143204      end
    144205
     206
     207     
     208      subroutine interpnir(escout,p,nlayer,escin,pin,nl)
     209C
     210C subroutine to perform linear interpolation in pressure from 1D profile
     211C escin(nl) sampled on pressure grid pin(nl) to profile
     212C escout(nlayer) on pressure grid p(nlayer).
     213C
     214      real escout(nlayer),p(nlayer)
     215      real escin(nl),pin(nl),wm,wp
     216      integer nl,nlayer,n1,n,nm,np
     217      do n1=1,nlayer
     218         if(p(n1) .gt. 1500. .or. p(n1) .lt. 1.0e-13) then
     219            escout(n1) = 0.0
     220         else
     221            do n = 1,nl-1
     222               if (p(n1).le.pin(n).and.p(n1).ge.pin(n+1)) then
     223                  nm=n
     224                  np=n+1
     225                  wm=abs(pin(np)-p(n1))/(pin(nm)-pin(np))
     226                  wp=1.0 - wm
     227               endif
     228            enddo
     229            escout(n1) = escin(nm)*wm + escin(np)*wp
     230         endif
     231      enddo
     232      return
     233      end

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

    r38 r414  
    11c**************************************************************************
    22c
    3       subroutine nltecool(ngrid,nlayer,pplay,pt,dtnlte)
     3      subroutine nltecool(ngrid,nlayer,nq,pplay,pt,pq,dtnlte)
    44c
    55c This code was designed as a delivery for the "Martian Environment Models"
    … …  
    2525c Version 1d  data contained in original input file "nlte_escape.dat"
    2626c now stored in include file "nltedata.h" Y.Wanherdrick 09/2000
     27
     28c       jul 2011 fgg   Modified to allow variable O
    2729c     
    2830c***************************************************************************
    … …  
    3133
    3234#include "nltedata.h" ! (Equivalent to the reading of the "nlte_escape.dat" file)
    33  
     35#include "dimensions.h"
     36#include "dimphys.h"
     37#include "chimiedata.h"
     38#include "conc.h" !Added to have "dynamic composition" in the scheme
     39#include "tracer.h" !"
     40#include "callkeys.h"
     41
    3442c Input and output variables
    3543c
    3644      integer     ngrid                           ! no. of horiz. gridpoints
    3745      integer     nlayer                          ! no. of atmospheric layers
     46      integer     nq                              ! no. of tracers
    3847      real        pplay(ngrid,nlayer)             ! input pressure grid
    3948      real        pt(ngrid,nlayer)                ! input temperatures
     49      real        pq(ngrid,nlayer,nq)                ! input mmrs
    4050      real        dtnlte(ngrid,nlayer)            ! output temp. tendencies
    4151
    … …  
    150160         call interp1(escf2,pyy,nlayer,ef2,pnb,np)
    151161         call interp1(escf1,pyy,nlayer,ef1,pnb,np)
    152          call interp3(co2,o3p,n2co,pyy,nlayer,
    153      &        co2vmr,o3pvmr,n2covmr,pnb,np)
     162         if(nltemodel.eq.0) then
     163            call interp3(co2,o3p,n2co,pyy,nlayer,
     164     &           co2vmr,o3pvmr,n2covmr,pnb,np)
     165         endif
    154166       
    155167         do i=1,nlayer  ! loop over layers
    … …  
    165177c           if(pt(j,i).lt.1.0)print*,pt(j,i)
    166178               nt = pyy(i)/(1.381e-17*pt(j,i)) ! nt in cm-3
     179               !Dynamic composition
     180               if(nltemodel.eq.1) then
     181                  co2(i)=pq(j,i,1)*mmean(j,i)/mmol(1)
     182                  o3p(i)=pq(j,i,3)*mmean(j,i)/mmol(3)
     183                  n2co(i)=pq(j,i,2)*mmean(j,i)/mmol(2) +
     184     $                 pq(j,i,12)*mmean(j,i)/mmol(12)
     185               endif
     186
     187               !Mixing ratio to density
    167188               co2(i)=co2(i)*nt                ! CO2 density in cm-3
    168189               o3p(i)=o3p(i)*nt                ! O3p density in cm-3

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

    r411 r414  
    6565c             Nb: See callradite.F for more information.
    6666c           Mesoscale lines: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags
     67c           jul 2011 malv+fgg: Modified calls to NIR heating routine and 15 um cooling parameterization
    6768c           
    6869c   arguments:
    … …  
    323324      REAL time_phys
    324325
     326! Added for new NLTE scheme
     327
     328      real co2vmr_gcm(ngridmx,nlayermx)
     329      real n2vmr_gcm(ngridmx,nlayermx)
     330      real ovmr_gcm(ngridmx,nlayermx)
     331      real covmr_gcm(ngridmx,nlayermx)
     332
     333
    325334c Variables for PBL
    326335
    … …  
    422431         endif         
    423432
     433         if(callnlte.and.nltemodel.eq.2) call NLTE_leedat
     434         if(callnirco2.and.nircorr.eq.1) call NIR_leedat
     435
    424436        IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN
    425437          write(*,*)"physiq: water_param Surface water ice albedo:",
    … …  
    428440                   
    429441      ENDIF        !  (end of "if firstcall")
     442
    430443
    431444c ---------------------------------------------------
    … …  
    532545c          NLTE cooling from CO2 emission
    533546c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    534 
    535            IF(callnlte) CALL nltecool(ngrid,nlayer,pplay,pt,zdtnlte)
     547           IF(callnlte) then
     548              if(nltemodel.eq.0.or.nltemodel.eq.1) then
     549                 CALL nltecool(ngrid,nlayer,nq,pplay,pt,pq,zdtnlte)
     550              else if(nltemodel.eq.2) then
     551                 do ig=1,ngrid
     552                    do l=1,nlayer
     553                       co2vmr_gcm(ig,l)=pq(ig,l,igcm_co2)*
     554     $                      mmean(ig,l)/mmol(igcm_co2)
     555                       n2vmr_gcm(ig,l)=pq(ig,l,igcm_n2)*
     556     $                      mmean(ig,l)/mmol(igcm_n2)
     557                       covmr_gcm(ig,l)=pq(ig,l,igcm_co)*
     558     $                      mmean(ig,l)/mmol(igcm_co)
     559                       ovmr_gcm(ig,l)=pq(ig,l,igcm_o)*
     560     $                      mmean(ig,l)/mmol(igcm_o)
     561                    enddo
     562                 enddo
     563                 
     564                 CALL NLTEdlvr09_TCOOL(ngrid,nlayer,pplay*9.869e-6,
     565     $                pt,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
     566     $                ovmr_gcm,  zdtnlte )
     567
     568                 do ig=1,ngrid
     569                    do l=1,nlayer
     570                       zdtnlte(ig,l)=zdtnlte(ig,l)/86400.
     571                    enddo
     572                 enddo
     573              endif
     574           endif
    536575
    537576c          Find number of layers for LTE radiation calculations
    … …  
    605644           zdtnirco2(:,:)=0
    606645           if (callnirco2) then
    607               call nirco2abs (ngrid,nlayer,pplay,dist_sol,
     646              call nirco2abs (ngrid,nlayer,pplay,dist_sol,nq,pq,
    608647     .                       mu0,fract,declin, zdtnirco2)
    609648           endif
    … …  
    9951034c        --------------
    9961035         IF (photochem .or. thermochem) then
     1036!NB: Photochemistry includes condensation of H2O2
    9971037            PRINT*, 'SURFDUST,SURFICE TO BE IMPLEMENTED. YAAAAAARG.'
    9981038            call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0,
Note: See TracChangeset for help on using the changeset viewer.