Changeset 414


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
Files:
54 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r412 r414  
    12441244
    12451245>>Commented aggressive outputs at first call in suaer.F90, aeropacity.F
     1246
     1247== 23/11/11 == FGG + MALV
     1248
     1249>> 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:
     1250
     1251-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.
     1252
     1253-inifis.F. Reading of nltemodel is added.
     1254
     1255-callkeys.h Declaration of nltemodel is added.
     1256
     1257The following lines should be added to callphys.def (ideally after setting callnlte):
     1258
     1259# NLTE 15um scheme to use.
     1260# 0->  Old scheme, static oxygen
     1261# 1->  Old scheme, dynamic oxygen
     1262# 2->  New scheme
     1263nltemodel = 2
     1264
     1265A new directory, NLTEDAT, has to be included in datagcm.
     1266
     1267>> 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:
     1268
     1269-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.
     1270
     1271-physiq.F: Call to NIR_leedat added at first call. Modified call to nirco2abs
     1272
     1273-inifis: Reading new flag nircorr.
     1274
     1275-callkeys.h: Declaration of nircorr.
     1276
     1277The following lines have to be added to callphys.def (ideally after callnirco2):
     1278
     1279# NIR NLTE correction for variable SZA and O/CO2?
     1280# matters only if callnirco2=T
     1281# 0-> no correction
     1282# 1-> include correction
     1283nircorr=1
     1284
     1285A new data file, NIRcorrection_feb2011.dat, has to be included in datagcm.
     1286
     1287>> 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

  • trunk/LMDZ.MARS/libf/aeronomars/moldiff.F

    r38 r414  
    3232c Local
    3333c
    34        real hco2(ncomptot),ho
     34
     35      integer,parameter :: ncompmoldiff = 14
     36       real hco2(ncompmoldiff),ho
    3537
    3638      integer ig,nz,l,n,nn,iq
    … …  
    3941      real hp(nlayermx)
    4042      real tt(nlayermx)
    41       real qq(nlayermx,ncomptot)
     43      real qq(nlayermx,ncompmoldiff)
    4244      real dmmeandz(nlayermx)
    43       real qnew(nlayermx,ncomptot)
     45      real qnew(nlayermx,ncompmoldiff)
    4446      real zlocal(nlayermx)
    45       real alf(ncomptot-1,ncomptot-1)
    46       real alfinv(nlayermx,ncomptot-1,ncomptot-1)
    47       real indx(ncomptot-1)
    48       real b(nlayermx,ncomptot-1)
    49       real y(ncomptot-1,ncomptot-1)
    50       real aa(nlayermx,ncomptot-1,ncomptot-1)
    51       real bb(nlayermx,ncomptot-1,ncomptot-1)
    52       real cc(nlayermx,ncomptot-1,ncomptot-1)
     47      real alf(ncompmoldiff-1,ncompmoldiff-1)
     48      real alfinv(nlayermx,ncompmoldiff-1,ncompmoldiff-1)
     49      real indx(ncompmoldiff-1)
     50      real b(nlayermx,ncompmoldiff-1)
     51      real y(ncompmoldiff-1,ncompmoldiff-1)
     52      real aa(nlayermx,ncompmoldiff-1,ncompmoldiff-1)
     53      real bb(nlayermx,ncompmoldiff-1,ncompmoldiff-1)
     54      real cc(nlayermx,ncompmoldiff-1,ncompmoldiff-1)
    5355      real atri(nlayermx-2)
    5456      real btri(nlayermx-2)
    … …  
    5658      real rtri(nlayermx-2)
    5759      real qtri(nlayermx-2)
    58       real alfdiag(ncomptot-1)
    59       real wi(ncomptot), flux(ncomptot), pote
     60      real alfdiag(ncompmoldiff-1)
     61      real wi(ncompmoldiff), flux(ncompmoldiff), pote
    6062
    6163cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    … …  
    9597      integer,save :: g_h2o=0
    9698
    97       integer,save :: gcmind(ncomptot) ! array of GCM indexes
     99      integer,save :: gcmind(ncompmoldiff) ! array of GCM indexes
    98100      integer ierr
    99101
    100102      logical,save :: firstcall=.true.
    101       real abfac(ncomptot)
    102       real,save :: dij(ncomptot,ncomptot)
     103      real abfac(ncompmoldiff)
     104      real,save :: dij(ncompmoldiff,ncompmoldiff)
    103105
    104106! Initializations at first call
    … …  
    215217     &                  +pdtconduc(ig,l)*ptimestep
    216218
    217           do nn=1,ncomptot
     219          do nn=1,ncompmoldiff
    218220            qq(l,nn)=pq(ig,l,gcmind(nn))+pdq(ig,l,gcmind(nn))*ptimestep
    219221            qq(l,nn)=max(qq(l,nn),1.e-30)
    … …  
    230232     &                  +pdtconduc(ig,nz)*ptimestep
    231233
    232         do nn=1,ncomptot
     234        do nn=1,ncompmoldiff
    233235          qq(1,nn)=pq(ig,1,gcmind(nn))+pdq(ig,1,gcmind(nn))*ptimestep
    234236          qq(nz,nn)=pq(ig,nz,gcmind(nn))+pdq(ig,nz,gcmind(nn))*ptimestep
    … …  
    251253        ntot=pplay(ig,l)/(1.381e-23*tt(l))                      ! in #/m3
    252254
    253         do nn=1,ncomptot-1            ! rows
     255        do nn=1,ncompmoldiff-1            ! rows
    254256          alfdiag(nn)=0.
    255257          dcoef1=dij(nn,i_o)*ptfac
    256           do n=1,ncomptot-1           ! columns
     258          do n=1,ncompmoldiff-1           ! columns
    257259            y(nn,n)=0.
    258260            dcoef=dij(nn,n)*ptfac
    … …  
    277279c Inverting the alfa matrix
    278280c
    279         call ludcmp(alf,ncomptot-1,ncomptot-1,indx,d,ierr)
     281        call ludcmp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,d,ierr)
    280282
    281283c       TEMPORAIRE *****************************
    … …  
    291293        end if
    292294c       *******************************************
    293         do n=1,ncomptot-1
    294           call lubksb(alf,ncomptot-1,ncomptot-1,indx,y(1,n))
    295           do nn=1,ncomptot-1
     295        do n=1,ncompmoldiff-1
     296          call lubksb(alf,ncompmoldiff-1,ncompmoldiff-1,indx,y(1,n))
     297          do nn=1,ncompmoldiff-1
    296298            alfinv(l,nn,n)=y(nn,n)/hh
    297299          enddo
    … …  
    309311        del1=hp(l)*pplay(ig,l)/(g*ptimestep)
    310312        del2=(hp(l)/2)**2*pplay(ig,l)/(g*ptimestep)
    311         do nn=1,ncomptot-1
    312           do n=1,ncomptot-1
     313        do nn=1,ncompmoldiff-1
     314          do n=1,ncompmoldiff-1
    313315            dalfinvdz=(alfinv(l+1,nn,n)-alfinv(l-1,nn,n))/hp(l)
    314316            aa(l,nn,n)=-dalfinvdz/del1+alfinv(l,nn,n)/del2
    … …  
    336338c (Hunten 1973, eq. 5)
    337339     
    338       do n=1,ncomptot
     340      do n=1,ncompmoldiff
    339341        wi(n)=1.
    340342        flux(n)=0.
    … …  
    379381c
    380382
    381       do nn=1,ncomptot-1
     383      do nn=1,ncompmoldiff-1
    382384        do l=2,nz-1
    383385          atri(l-1)=aa(l,nn,nn)
    … …  
    385387          ctri(l-1)=cc(l,nn,nn)
    386388          rtri(l-1)=qq(l,nn)
    387           do n=1,ncomptot-1
     389          do n=1,ncompmoldiff-1
    388390            rtri(l-1)=rtri(l-1)-(aa(l,nn,n)*qq(l-1,n)
    389391     &                           +bb(l,nn,n)*qq(l,n)
    … …  
    434436        if(zlocal(l).gt.65000.)then
    435437        pdqdiff(ig,l,g_o)=0.
    436         do n=1,ncomptot-1
     438        do n=1,ncompmoldiff-1
    437439          pdqdiff(ig,l,gcmind(n))=(qnew(l,n)-qq(l,n))
    438440          pdqdiff(ig,l,g_o)=pdqdiff(ig,l,g_o)-(qnew(l,n)-qq(l,n))

  • trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff.F

    r38 r414  
    2222c    Input/Output
    2323c    ------------
    24       real dij(ncomptot,ncomptot)
     24       integer,parameter :: ncompmoldiff = 14
     25      real dij(ncompmoldiff,ncompmoldiff)
    2526
    2627c    Local variables:
    … …  
    6364      integer,save :: g_h2o=0
    6465
    65       integer,save :: gcmind(ncomptot)
     66      integer,save :: gcmind(ncompmoldiff)
    6667
    6768      real dnh
    … …  
    187188       print*,'moldiffcoef: COEFF CALC'
    188189       open(56,file='coeffs.dat',status='unknown')
    189       do n=1,ncomptot
     190      do n=1,ncompmoldiff
    190191        if (dij(i_h2,n).gt.0.0) then
    191           do nn=n,ncomptot
     192          do nn=n,ncompmoldiff
    192193            dij(nn,n)=dij(i_h2,n)
    193194     &                  *sqrt(mmol(g_h2)/mmol(gcmind(nn)))
    … …  
    198199        if (dij(i_h2,n).eq.0.0) then
    199200          dnh=dij(i_h,i_o)*sqrt(mmol(g_o)/mmol(gcmind(n)))
    200           do nn=n,ncomptot
     201          do nn=n,ncompmoldiff
    201202            dij(nn,n)=dnh*sqrt(mmol(g_h)/mmol(gcmind(nn)))
    202203            if(n.eq.nn) dij(nn,n)=1.0
    … …  
    206207      enddo
    207208
    208       do n=1,ncomptot
    209         do nn=n,ncomptot
     209      do n=1,ncompmoldiff
     210        do nn=n,ncompmoldiff
    210211          write(56,*) n,nn,dij(n,nn)    !*1.e5/1.381e-23/(273**1.75)
    211212        enddo

  • 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.