Changeset 2464 for trunk


Ignore:
Timestamp:
Feb 22, 2021, 4:43:11 PM (4 years ago)
Author:
slebonnois
Message:

SL: major update related to extension to 250 km. New EUV heating as used in Martian GCM, debug of a few routines, adding He, clean up, updating mmean regularly, modification of the order of processes, add the possibility to use Hedin composition in rad transfer

Location:
trunk/LMDZ.VENUS/libf/phyvenus
Files:
1 added
4 deleted
20 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.VENUS/libf/phyvenus/chemparam_mod.F90

    r2048 r2464  
    1313                 i_cocl2, i_s, i_so, i_so2, i_so3,       &
    1414                 i_s2o2, i_ocs, i_hso3, i_h2so4, i_s2,   &
    15                  i_clso2, i_oscl, i_n2
     15                 i_clso2, i_oscl, i_n2, i_he
    1616                 
    1717INTEGER, SAVE :: i_h2oliq, i_h2so4liq
    … …  
    834834                i_n2=i
    835835                M_tr(i_n2)=28.013
     836                CASE('he')
     837                i_he=i
     838                M_tr(i_he)=4.0026
    836839! MICROPHYSICAL TRACERS FOR CL_SCHEME=1
    837840                CASE('h2oliq')

  • trunk/LMDZ.VENUS/libf/phyvenus/clesphys.h

    r2135 r2464  
    2323       REAL    z0, lmixmin
    2424       REAL    ksta, inertie
    25        REAL    euveff, solarcondate
     25       REAL    euveff, fixed_euv_value
    2626
    2727       COMMON/clesphys_l/ cycle_diurne, soil_model,                     &
    … …  
    3737
    3838       COMMON/clesphys_r/ ecriphy, solaire, z0, lmixmin,                &
    39      &     ksta, inertie, euveff, solarcondate
     39     &     ksta, inertie, euveff, fixed_euv_value
    4040

  • trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/new_cloud_sedim.F

    r2203 r2464  
    9696     $            + d_tr_chem(:,:,i_h2so4liq)*ptimestep
    9797
    98       wgt_SA(:,:) = wh2so4(:,:)
     98      wgt_SA(:,1:n_lev) = wh2so4(:,:)
    9999
    100100c     Init F_sed

  • trunk/LMDZ.VENUS/libf/phyvenus/compo_hedin_mod2.F90

    r1310 r2464  
    9797         enddo
    9898
    99 !            print*, pres_hedin(1,:)
    100 !           print*, ' '
     99!           print*, pres_hedin(1,:)
     100!           print*, ' '
    101101!           print*, T(musize,:) 
    102102!           print*, ' '
    … …  
    112112!               print*, " "
    113113!               print*, ntot(1,:)
    114 !               print*, co2_hedin(:,10)
     114!               print*, co2_hedin(1,:)
    115115!               print*, mu_hedin(:)
    116 !               print*, " "
    117 !               print*, n2_hedin(1,:)
    118 !               print*, " "
    119 !               print*, o_hedin(1,:)
    120 !               print*, " "
    121 !               print*, co_hedin(1,:)
    122 !               print*, ' '
    123 !        stop
    124                
     116!               print*, " "
     117!               print*, n2_hedin(1,:)
     118!               print*, " "
     119!               print*, o_hedin(1,:)
     120!               print*, " "
     121!               print*, co_hedin(1,:)
     122!               print*, ' '
     123!          stop
     124
    125125       end subroutine compo_hedin83_init2
    126126
    … …  
    165165          enddo
    166166
    167         factp(k) = (log10(pression(i,k))-log10(pres_hedin(mu_ok,z_ok-1)))/(log10(pres_hedin(mu_ok,z_ok))-log10(pres_hedin(mu_ok,z_ok-1)))
     167        if (pres_hedin(mu_ok,zsize).ge.pression(i,k)) then
     168! avoid the extrapolation... Problem when p < 1.E-6 Pa...
     169         factp(k) = 1.
     170        else
     171         factp(k) = (log10(pression(i,k))-log10(pres_hedin(mu_ok,z_ok-1)))/(log10(pres_hedin(mu_ok,z_ok))-log10(pres_hedin(mu_ok,z_ok-1)))
     172        endif
    168173
    169        
    170174        ovmr_gcm(i,k) = (10.**(log10(o_hedin(mu_ok,z_ok))*factp(k)+log10(o_hedin(mu_ok,z_ok-1))*(1.-factp(k))))
    171175        n2vmr_gcm(i,k) = 10.**(log10(n2_hedin(mu_ok,z_ok))*factp(k)+log10(n2_hedin(mu_ok,z_ok-1))*(1.-factp(k)))

  • trunk/LMDZ.VENUS/libf/phyvenus/conc.F90

    r1525 r2464  
    1414      REAL, ALLOCATABLE, SAVE :: cpnew(:,:)
    1515!$OMP THREADPRIVATE(cpnew)
    16       REAL, ALLOCATABLE, SAVE :: jfotsout(:,:,:)
    17 !$OMP THREADPRIVATE(jfotsout)
    18       REAL, ALLOCATABLE, SAVE :: coefit4(:,:)
    19 !$OMP THREADPRIVATE(coefit4)
    20       REAL, ALLOCATABLE, SAVE :: coefit3(:,:)
    21 !$OMP THREADPRIVATE(coefit3)
    22       REAL, ALLOCATABLE, SAVE :: coefit2(:,:)
    23 !$OMP THREADPRIVATE(coefit2)
    24       REAL, ALLOCATABLE, SAVE :: coefit1(:,:)
    25 !$OMP THREADPRIVATE(coefit1)
    26       REAL, ALLOCATABLE, SAVE :: coefit0(:,:)
    27 !$OMP THREADPRIVATE(coefit0)
    28       REAL, ALLOCATABLE, SAVE :: fluxtop(:)
    29 !$OMP THREADPRIVATE(fluxtop)
    30 !      REAL, ALLOCATABLE, SAVE :: freccen(:)
    31 !$OMP THREADPRIVATE(freccen)
    32       REAL, ALLOCATABLE, SAVE :: t0(:)
    33 !$OMP THREADPRIVATE(t0)
    34 !      REAL, ALLOCATABLE, SAVE :: crscabsi2(:,:)
    35 !$OMP THREADPRIVATE(crscabsi2)
    36       REAL, ALLOCATABLE, SAVE :: jabsifotsintpar(:,:,:)
    37 !$OMP THREADPRIVATE(crscabsi2)
    38       REAL, ALLOCATABLE, SAVE :: c1_16(:,:)
    39 !$OMP THREADPRIVATE(c1_16)
    40       REAL, ALLOCATABLE, SAVE :: c17_24(:)
    41 !$OMP THREADPRIVATE(c17_24)
    42       REAL, ALLOCATABLE, SAVE :: c25_29(:)
    43 !$OMP THREADPRIVATE(c25_29)
    44       REAL, ALLOCATABLE, SAVE :: c30_31(:)
    45 !$OMP THREADPRIVATE(c30_31)
    46       REAL, ALLOCATABLE, SAVE :: c32(:)
    47 !$OMP THREADPRIVATE(c32)
    48       REAL, ALLOCATABLE, SAVE :: c33(:)
    49 !$OMP THREADPRIVATE(c33)
    50       REAL, ALLOCATABLE, SAVE :: c34(:)
    51 !$OMP THREADPRIVATE(c34)
    52       REAL, ALLOCATABLE, SAVE :: c35(:)
    53 !$OMP THREADPRIVATE(c35)
    54       REAL, ALLOCATABLE, SAVE :: c36(:)
    55 !$OMP THREADPRIVATE(c36)
    56       REAL, ALLOCATABLE, SAVE :: ct1(:)
    57 !$OMP THREADPRIVATE(ct1)
    58       REAL, ALLOCATABLE, SAVE :: ct2(:)
    59 !$OMP THREADPRIVATE(ct2)
    60       REAL, ALLOCATABLE, SAVE :: p1(:)
    61 !$OMP THREADPRIVATE(p1)
    62       REAL, ALLOCATABLE, SAVE :: p2(:)
    63 !$OMP THREADPRIVATE(p2)
    6416
    6517
    … …  
    6921     USE dimphy
    7022     implicit none
    71 #include"param.h"
    7223
    7324     ALLOCATE(mmean(klon,klev))
    … …  
    7627     ALLOCATE(rnew(klon,klev)) 
    7728     ALLOCATE(cpnew(klon,klev))
    78      ALLOCATE(jfotsout(ninter,nabs,klev))
    79      ALLOCATE(coefit0(ninter,nabs))
    80      ALLOCATE(coefit1(ninter,nabs)) 
    81      ALLOCATE(coefit2(ninter,nabs)) 
    82      ALLOCATE(coefit3(ninter,nabs)) 
    83      ALLOCATE(coefit4(ninter,nabs))
    84      ALLOCATE(fluxtop(ninter))
    85 !     ALLOCATE(freccen(ninter))
    86      ALLOCATE(t0(nz2))   
    87 !     ALLOCATE(crscabsi2(nabs,16))
    88      ALLOCATE(c1_16(nz2,16))
    89      ALLOCATE(c17_24(nz2) )
    90      ALLOCATE(c25_29(nz2) )
    91      ALLOCATE(c30_31(nz2) )
    92      ALLOCATE(c32(nz2) )
    93      ALLOCATE(c33(nz2) )
    94      ALLOCATE(c34(nz2) )
    95      ALLOCATE(c35(nz2) )
    96      ALLOCATE(c36(nz2) )
    97      ALLOCATE(ct2(ninter) )
    98      ALLOCATE(ct1(ninter) )
    99      ALLOCATE(p1(ninter) )
    100      ALLOCATE(p2(ninter) )
    101      ALLOCATE(jabsifotsintpar(nz2,nabs,ninter))
    10229     end subroutine conc_init
    10330

  • trunk/LMDZ.VENUS/libf/phyvenus/concentrations2.F

    r1621 r2464  
    1       SUBROUTINE concentrations2(pplay,t_seri,pdt,tr_seri, nqmx)
     1      SUBROUTINE concentrations2(pplay,t_seri,tr_seri, nqmx)
    22
    33      use dimphy
    … …  
    3131
    3232      real pplay(klon,klev)
    33 c      real pt(klon,klev)
    3433      integer,intent(in) :: nqmx    ! number of tracers
    3534      real t_seri(klon, klev)
    36       real pdt(klon,klev)
    37       real n2vmr_gcm(klon,klev),nvmr_gcm(klon,klev)
    3835      real tr_seri(klon,klev,nqmx)
    39 c      real pdq(klon,klev,nqmx)
    40       real ptimestep
    4136
    4237!     local variables

  • trunk/LMDZ.VENUS/libf/phyvenus/conf_phys.F90

    r2135 r2464  
    454454
    455455!
    456 !Config Key  = solarcondate
     456!Config Key  = fixed_euv_value
    457457!Config Desc =
    458 !Config Def  = 1993.4  ## Average solar cycle condition
    459 !Config Help =
    460 !
    461   solarcondate = 1993.4
    462   call getin('solarcondate',solarcondate)
     458!Config Def  = 140.  ## Average solar cycle condition
     459!Config Help =
     460!
     461  fixed_euv_value =140.
     462  call getin('fixed_euv_value',fixed_euv_value)
    463463
    464464!
    … …  
    528528  write(lunout,*)' callthermos = ',callthermos
    529529  write(lunout,*)' solvarmod = ',solvarmod
    530   write(lunout,*)' solarcondate = ',solarcondate
     530  write(lunout,*)' fixed_euv_value = ',fixed_euv_value
    531531  write(lunout,*)' euveff = ',euveff
    532532

  • trunk/LMDZ.VENUS/libf/phyvenus/euvheat.F90

    r1530 r2464  
    11      SUBROUTINE euvheat(nlon, nlev,nqmx, pt,pplev,pplay,zzlay, &
    2                    mu0,ptimestep,ptime,zday,       &
    3                    pq, pdq, pdteuv)
     2                   mu0,ptimestep,ptime,pq, pdq, pdteuv)
    43
    54        use chemparam_mod
    6         use dimphy
     5        use dimphy
    76        use conc, only:  rnew, cpnew
    87
    … …  
    3433#include "YOMCST.h"
    3534#include "clesphys.h"
    36 #include "param.h"
    37 #include "param_v4.h"
    38 !#include "chimiedata.h"
    3935#include "mmol.h"
    4036!-----------------------------------------------------------------------
    … …  
    5450
    5551      real :: ptimestep,ptime
    56       real :: zday
    5752      real :: pq(nlon,nlev,nqmx)
    5853      real :: pdq(nlon,nlev,nqmx)
    … …  
    336331     
    337332      !Solar flux calculation
    338       if(solvarmod.eq.0) call flujo(solarcondate)
    339 
    340                                          ! Not recommended for long runs
    341                                          ! (e.g. to build the MCD), as the
    342                                          ! solar conditions at the end will
    343                                          ! be different to the ones initially
    344                                          ! set
    345333     
    346334      do ig=1,nlon
    … …  
    377365
    378366        !Routine to calculate the UV heating
    379          call hrtherm (ig,euvmod,rm,nespeuv,tx,zlocal,zenit,zday,jtot)
     367         call hrtherm (ig,euvmod,rm,nespeuv,tx,zlocal,zenit,jtot)
    380368         
    381369

  • trunk/LMDZ.VENUS/libf/phyvenus/hrtherm.F

    r1530 r2464  
    11c**********************************************************************
    22
    3       subroutine hrtherm(ig,euvmod,rm,nespeuv,tx,iz,zenit,zday,jtot)
     3      subroutine hrtherm(ig,euvmod,rm,nespeuv,tx,iz,zenit,jtot)
    44
    55
    … …  
    99c**********************************************************************
    1010      use dimphy
    11       use conc
     11      use param_v4_h, only: ninter,nabs,jfotsout,fluxtop,freccen
     12
    1213      implicit none
    1314
    … …  
    1516
    1617
    17 #include "param.h"
    18 #include "param_v4.h"
    1918#include "clesphys.h"
    2019
    … …  
    3837      real       zenit
    3938      real       iz(klev)
    40       real       zday
    4139
    4240      ! tracer indexes for the EUV heating:
    … …  
    9593
    9694      !Calculation of photoabsortion coefficient
    97       if(solvarmod.eq.0) then
    98          call jthermcalc(ig,euvmod,rm,nespeuv,tx,iz,zenit)
    99       else if (solvarmod.eq.1) then
    100          call jthermcalc_e107(ig,euvmod,rm,nespeuv,tx,iz,zenit,zday)
    101          do indexint=1,ninter
    102             fluxtop(indexint)=1.
    103          enddo
    104       endif
     95      call jthermcalc_e107(ig,klev,euvmod,rm,nespeuv,tx,iz,zenit)
    10596
    10697      !Total photoabsorption coefficient    !  erg/(s*cm3)

  • trunk/LMDZ.VENUS/libf/phyvenus/jthermcalc_e107.F

    r1530 r2464  
    22
    33      subroutine jthermcalc_e107
    4      $     (ig,chemthermod,rm,nesptherm,tx,iz,zenit,zday)
     4     $     (ig,nlayer,chemthermod,rm,nesptherm,tx,iz,zenit)
    55
    66
    … …  
    1111c MAC July 2003
    1212c**********************************************************************
    13         use dimphy
    14         use conc
     13
     14      use param_v4_h, only: jfotsout,crscabsi2,
     15     .    c1_16,c17_24,c25_29,c30_31,c32,c33,c34,c35,c36,
     16     .    co2crsc195,co2crsc295,t0,
     17     .    jabsifotsintpar,ninter,nz2,
     18     .    nabs,e107,date_e107,e107_tab,
     19     .    coefit0,coefit1,coefit2,coefit3,coefit4
     20
    1521      implicit none
    1622
    17 c     common variables and constants
    18 #include "param.h"
    19 #include "param_v4.h"
     23      include "clesphys.h"
     24
    2025
    2126c     input and output variables
    2227
    23       integer    ig
    24       integer    chemthermod
    25       integer    nesptherm                      !Number of species considered
    26       real       rm(klev,nesptherm)         !Densities (cm-3)
    27       real       tx(klev)                   !temperature
    28       real       zenit                          !SZA
    29       real       iz(klev)                   !Local altitude
    30       real       zday                           !Martian day after Ls=0
    31 
     28      integer,intent(in) :: ig,nlayer
     29      integer,intent(in) :: chemthermod
     30      integer,intent(in) :: nesptherm          !Number of species considered
     31      real,intent(in) :: rm(nlayer,nesptherm)  !Densities (cm-3)
     32      real,intent(in) :: tx(nlayer)            !temperature
     33      real,intent(in) :: zenit         !SZA
     34      real,intent(in) :: iz(nlayer)    !Local altitude
     35
     36! NB: no output variable! (computed jfotsout() is in module param_v4_h)
    3237
    3338c    local parameters and variables
    3439
    35       real       co2colx(klev)              !column density of CO2 (cm^-2)
    36       real       o2colx(klev)               !column density of O2(cm^-2)
    37       real       o3pcolx(klev)              !column density of O(3P)(cm^-2)
    38       real       h2colx(klev)               !H2 column density (cm-2)
    39       real       h2ocolx(klev)              !H2O column density (cm-2)
    40       real       h2o2colx(klev)             !column density of H2O2(cm^-2)
    41       real       o3colx(klev)               !O3 column density (cm-2)
    42       real       n2colx(klev)               !N2 column density (cm-2)
    43       real       ncolx(klev)                !N column density (cm-2)
    44       real       nocolx(klev)               !NO column density (cm-2)
    45       real       cocolx(klev)               !CO column density (cm-2)
    46       real       hcolx(klev)                !H column density (cm-2)
    47       real       no2colx(klev)              !NO2 column density (cm-2)
    48       real       t2(klev)
    49       real       coltemp(klev)
    50       real       sigma(ninter,klev)
    51       real       alfa(ninter,klev)
    52       real       realday
     40      real, parameter :: dist_sol=0.72333
     41
     42      real       co2colx(nlayer)              !column density of CO2 (cm^-2)
     43      real       o2colx(nlayer)               !column density of O2(cm^-2)
     44      real       o3pcolx(nlayer)              !column density of O(3P)(cm^-2)
     45      real       h2colx(nlayer)               !H2 column density (cm-2)
     46      real       h2ocolx(nlayer)              !H2O column density (cm-2)
     47      real       h2o2colx(nlayer)             !column density of H2O2(cm^-2)
     48      real       o3colx(nlayer)               !O3 column density (cm-2)
     49      real       n2colx(nlayer)               !N2 column density (cm-2)
     50      real       ncolx(nlayer)                !N column density (cm-2)
     51      real       nocolx(nlayer)               !NO column density (cm-2)
     52      real       cocolx(nlayer)               !CO column density (cm-2)
     53      real       hcolx(nlayer)                !H column density (cm-2)
     54      real       no2colx(nlayer)              !NO2 column density (cm-2)
     55      real       t2(nlayer)
     56      real       coltemp(nlayer)
     57      real       sigma(ninter,nlayer)
     58      real       alfa(ninter,nlayer)
    5359     
    5460      integer    i,j,k,indexint                 !indexes
    … …  
    7480      real*8      auxjh(nz2)
    7581      real*8      auxjno2(nz2)
    76       real*8      wp(klev),wm(klev)
    77       real*8      auxcolinp(klev)
    78       integer     auxind(klev)
     82      real*8      wp(nlayer),wm(nlayer)
     83      real*8      auxcolinp(nlayer)
     84      integer     auxind(nlayer)
    7985      integer     auxi
    8086      integer     ind
    81       real*8      cortemp(klev)
     87      real*8      cortemp(nlayer)
    8288
    8389      real*8      limdown                      !limits for interpolation
    8490      real*8      limup                        !        ""
    8591
    86       !!!ATTENTION. Here i_co2 has to have the same value than in chemthermos.F90
     92      !!!ATTENTION. Here i_co2 has to have the same value than in euvheat.F90
    8793      !!!If the value is changed there, if has to be changed also here !!!!
    8894      integer,parameter :: i_co2=1
    8995
     96      character*20 modname
     97      character*80 abort_message
    9098
    9199c*************************PROGRAM STARTS*******************************
    92100     
     101      modname = 'jthermcalc_e107'
     102
    93103      if(zenit.gt.140.) then
    94104         dn='n'
    … …  
    105115      !Auxiliar temperature to take into account the temperature dependence
    106116      !of CO2 cross section
    107       do i=1,klev
     117      do i=1,nlayer
    108118         t2(i)=tx(i)
    109119         if(t2(i).lt.195.0) t2(i)=195.0
    … …  
    113123      !Calculation of column amounts
    114124      call column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
    115      $     co2colx,o2colx,o3pcolx,h2colx,h2ocolx,
    116      $     h2o2colx,o3colx,n2colx,ncolx,nocolx,cocolx,hcolx,no2colx)
     125     $     co2colx,o3pcolx,n2colx,cocolx)
    117126
    118127      !Auxiliar column to include the temperature dependence
    119128      !of CO2 cross section
    120       coltemp(klev)=co2colx(klev)*abs(t2(klev)-t0(klev))
    121       do i=klev-1,1,-1
     129      coltemp(nlayer)=co2colx(nlayer)*abs(t2(nlayer)-t0(nlayer))
     130      do i=nlayer-1,1,-1
    122131        coltemp(i)=!coltemp(i+1)+     PQ SE ELIMINA? REVISAR
    123132     $         ( rm(i,i_co2) + rm(i+1,i_co2) ) * 0.5
    … …  
    126135     
    127136      !Calculation of CO2 cross section at temperature t0(i)
    128       do i=1,klev
     137      do i=1,nlayer
    129138         do indexint=24,32
    130139           sigma(indexint,i)=co2crsc195(indexint-23)
    … …  
    134143      end do
    135144
    136       !E10.7 for the day: linear interpolation to tabulated values
    137       realday=mod(zday,669.)
    138       if(realday.lt.date_e107(1)) then
    139          e107=e107_tab(1)
    140       else if(realday.ge.date_e107(669)) then
    141          e107=e107_tab(669)   
    142       else if(realday.ge.date_e107(1).and.
    143      $        realday.lt.date_e107(669)) then
    144          do i=1,668
    145             if(realday.ge.date_e107(i).and.
    146      $           realday.lt.date_e107(i+1)) then
    147                tinf=i
    148                tsup=i+1
    149                e107=e107_tab(tinf)+(realday-date_e107(tinf))*
    150      $              (e107_tab(tsup)-e107_tab(tinf))
    151             endif
    152          enddo
    153       endif
     145      if (solvarmod==0) then
     146        e107=fixed_euv_value
     147      else
     148            abort_message='solvarmod should be 0...'
     149            call abort_physic(modname,abort_message,1)
     150      endif ! of if (solvarmod==0)
    154151
    155152      !Photoabsorption coefficients at TOA as a function of E10.7
    156153      do j=1,nabs
    157154         do indexint=1,ninter
    158             jfotsout(indexint,j,klev)=coefit0(indexint,j)+
     155            jfotsout(indexint,j,nlayer)=coefit0(indexint,j)+
    159156     $           coefit1(indexint,j)*e107+coefit2(indexint,j)*e107**2+
    160157     $           coefit3(indexint,j)*e107**3+coefit4(indexint,j)*e107**4
    … …  
    178175c     Input atmospheric column
    179176      indexint=1
    180       do i=1,klev
    181          auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint) +
     177      do i=1,nlayer
     178         auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint) +
    182179     $        o2colx(i)*crscabsi2(2,indexint) +
    183180     $        o3pcolx(i)*crscabsi2(3,indexint) +
    … …  
    213210
    214211      call interfast
    215      $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
    216       do i=1,klev
     212     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
     213      do i=1,nlayer
    217214         ind=auxind(i)
    218          auxi=klev-i+1
     215         auxi=nlayer-i+1
    219216         !CO2 interpolated coefficient
    220          jfotsout(indexint,1,auxi) = jfotsout(indexint,1,klev) *
     217         jfotsout(indexint,1,auxi) = jfotsout(indexint,1,nlayer) *
    221218     $        (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
    222219         !O2 interpolated coefficient
    223          jfotsout(indexint,2,auxi) = jfotsout(indexint,2,klev) *
     220         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
    224221     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
    225222         !O3p interpolated coefficient
    226          jfotsout(indexint,3,auxi) = jfotsout(indexint,3,klev) *
     223         jfotsout(indexint,3,auxi) = jfotsout(indexint,3,nlayer) *
    227224     $        (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
    228225         !H2 interpolated coefficient
    … …  
    233230      !N interpolated coefficient
    234231      if(chemthermod.ge.2) then
    235          do i=1,klev
     232         do i=1,nlayer
    236233            ind=auxind(i)
    237             jfotsout(indexint,9,klev-i+1) = 
    238      $           jfotsout(indexint,9,klev) *
     234            jfotsout(indexint,9,nlayer-i+1) = 
     235     $           jfotsout(indexint,9,nlayer) *
    239236     $           (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
    240237         enddo
    … …  
    255252c     Input atmospheric column
    256253      do indexint=2,15
    257          do i=1,klev
    258             auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint)+
     254         do i=1,nlayer
     255            auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint)+
    259256     $           o2colx(i)*crscabsi2(2,indexint)+
    260257     $           o3pcolx(i)*crscabsi2(3,indexint)+
    … …  
    302299         endif
    303300
    304           call interfast(wm,wp,auxind,auxcolinp,klev,
     301          call interfast(wm,wp,auxind,auxcolinp,nlayer,
    305302     $        auxcoltab,nz2,limdown,limup)
    306           do i=1,klev
     303          do i=1,nlayer
    307304             ind=auxind(i)
    308              auxi = klev-i+1
     305             auxi = nlayer-i+1
    309306             !O2 interpolated coefficient
    310307             jfotsout(indexint,2,auxi) =
    311      $            jfotsout(indexint,2,klev) *
     308     $            jfotsout(indexint,2,nlayer) *
    312309     $            (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
    313310             !O3p interpolated coefficient
    314311             jfotsout(indexint,3,auxi) =
    315      $            jfotsout(indexint,3,klev) *
     312     $            jfotsout(indexint,3,nlayer) *
    316313     $            (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
    317314             !CO2 interpolated coefficient
    318315             jfotsout(indexint,1,auxi) =
    319      $            jfotsout(indexint,1,klev) *
     316     $            jfotsout(indexint,1,nlayer) *
    320317     $            (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
    321318             !H2 interpolated coefficient
    322319             jfotsout(indexint,5,auxi) =
    323      $            jfotsout(indexint,5,klev) *
     320     $            jfotsout(indexint,5,nlayer) *
    324321     $            (wm(i)*auxjh2(ind+1) + wp(i)*auxjh2(ind))
    325322             !N2 interpolated coefficient
    326323             jfotsout(indexint,8,auxi) =
    327      $            jfotsout(indexint,8,klev) *
     324     $            jfotsout(indexint,8,nlayer) *
    328325     $            (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind))
    329326             !CO interpolated coefficient
    330327             jfotsout(indexint,11,auxi) =
    331      $            jfotsout(indexint,11,klev) *
     328     $            jfotsout(indexint,11,nlayer) *
    332329     $            (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind))
    333330             !H interpolated coefficient
    334331             jfotsout(indexint,12,auxi) =
    335      $            jfotsout(indexint,12,klev) *
     332     $            jfotsout(indexint,12,nlayer) *
    336333     $            (wm(i)*auxjh(ind+1) + wp(i)*auxjh(ind))
    337334          enddo
    338335          !Only if chemthermod.ge.2
    339336          if(chemthermod.ge.2) then
    340              do i=1,klev
     337             do i=1,nlayer
    341338                ind=auxind(i)
    342                 auxi = klev-i+1
     339                auxi = nlayer-i+1
    343340                !N interpolated coefficient
    344341                jfotsout(indexint,9,auxi) =
    345      $               jfotsout(indexint,9,klev) *
     342     $               jfotsout(indexint,9,nlayer) *
    346343     $               (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
    347344                !NO interpolated coefficient
    348345                jfotsout(indexint,10,auxi)=
    349      $               jfotsout(indexint,10,klev) *
     346     $               jfotsout(indexint,10,nlayer) *
    350347     $               (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind))
    351348                !NO2 interpolated coefficient
    352349                jfotsout(indexint,13,auxi)=
    353      $               jfotsout(indexint,13,klev) *
     350     $               jfotsout(indexint,13,nlayer) *
    354351     $               (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
    355352             enddo
    … …  
    370367c     Input atmospheric column
    371368      indexint=16
    372       do i=1,klev
    373          auxcolinp(klev-i+1) = co2colx(i)*crscabsi2(1,indexint)+
     369      do i=1,nlayer
     370         auxcolinp(nlayer-i+1) = co2colx(i)*crscabsi2(1,indexint)+
    374371     $        o2colx(i)*crscabsi2(2,indexint)+
    375372     $        o3pcolx(i)*crscabsi2(3,indexint)+
    … …  
    417414
    418415      call interfast
    419      $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
    420       do i=1,klev
     416     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
     417      do i=1,nlayer
    421418         ind=auxind(i)
    422          auxi = klev-i+1
     419         auxi = nlayer-i+1
    423420         !O2 interpolated coefficient
    424          jfotsout(indexint,2,auxi) = jfotsout(indexint,2,klev) *
     421         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
    425422     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
    426423         !CO2 interpolated coefficient
    427          jfotsout(indexint,1,auxi) = jfotsout(indexint,1,klev) *
     424         jfotsout(indexint,1,auxi) = jfotsout(indexint,1,nlayer) *
    428425     $        (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
    429426         !O3p interpolated coefficient
    430          jfotsout(indexint,3,auxi) = jfotsout(indexint,3,klev) *
     427         jfotsout(indexint,3,auxi) = jfotsout(indexint,3,nlayer) *
    431428     $        (wm(i)*auxjo3p(ind+1) + wp(i)*auxjo3p(ind))
    432429         !N2 interpolated coefficient
    433          jfotsout(indexint,8,auxi) = jfotsout(indexint,8,klev) *
     430         jfotsout(indexint,8,auxi) = jfotsout(indexint,8,nlayer) *
    434431     $        (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind))
    435432         !CO interpolated coefficient
    436433         jfotsout(indexint,11,auxi) =
    437      $        jfotsout(indexint,11,klev) *
     434     $        jfotsout(indexint,11,nlayer) *
    438435     $        (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind))
    439436         !H interpolated coefficient
    440437         jfotsout(indexint,12,auxi) =
    441      $        jfotsout(indexint,12,klev) *
     438     $        jfotsout(indexint,12,nlayer) *
    442439     $        (wm(i)*auxjh(ind+1) + wp(i)*auxjh(ind))
    443440      enddo
    444441      !Only if chemthermod.ge.2
    445442      if(chemthermod.ge.2) then
    446          do i=1,klev
     443         do i=1,nlayer
    447444            ind=auxind(i)
    448             auxi = klev-i+1
     445            auxi = nlayer-i+1
    449446            !N interpolated coefficient
    450447            jfotsout(indexint,9,auxi) =
    451      $           jfotsout(indexint,9,klev) *
     448     $           jfotsout(indexint,9,nlayer) *
    452449     $           (wm(i)*auxjn(ind+1) + wp(i)*auxjn(ind))
    453450            !NO interpolated coefficient
    454451            jfotsout(indexint,10,auxi) =
    455      $           jfotsout(indexint,10,klev) *
     452     $           jfotsout(indexint,10,nlayer) *
    456453     $           (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind))
    457454            !NO2 interpolated coefficient
    458455            jfotsout(indexint,13,auxi) =
    459      $           jfotsout(indexint,13,klev) *
     456     $           jfotsout(indexint,13,nlayer) *
    460457     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
    461458         enddo
    … …  
    473470c     Input column
    474471
    475       do i=1,klev
    476          auxcolinp(klev-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
     472      do i=1,nlayer
     473         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + n2colx(i) +
    477474     $        nocolx(i) + cocolx(i) + no2colx(i)
    478475      end do
    … …  
    507504
    508505         call interfast
    509      $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
     506     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
    510507         !Correction to include T variation of CO2 cross section
    511508         if(indexint.eq.24) then
    512             do i=1,klev
    513                auxi = klev-i+1
     509            do i=1,nlayer
     510               auxi = nlayer-i+1
    514511               if(sigma(indexint,auxi)*
    515512     $              alfa(indexint,auxi)*coltemp(auxi)
    … …  
    522519            enddo
    523520         else
    524             do i=1,klev
     521            do i=1,nlayer
    525522               cortemp(i)=1.
    526523            enddo
    527524         end if
    528          do i=1,klev           
     525         do i=1,nlayer           
    529526            ind=auxind(i)
    530             auxi = klev-i+1
     527            auxi = nlayer-i+1
    531528            !O2 interpolated coefficient
    532529            jfotsout(indexint,2,auxi) =
    533      $           jfotsout(indexint,2,klev) *
     530     $           jfotsout(indexint,2,nlayer) *
    534531     $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) *
    535532     $           cortemp(i)
    536533            !CO2 interpolated coefficient
    537534            jfotsout(indexint,1,auxi) =
    538      $           jfotsout(indexint,1,klev) *
     535     $           jfotsout(indexint,1,nlayer) *
    539536     $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind))
    540537     $           * cortemp(i)
    … …  
    545542            !N2 interpolated coefficient
    546543            jfotsout(indexint,8,auxi) =
    547      $           jfotsout(indexint,8,klev) *
     544     $           jfotsout(indexint,8,nlayer) *
    548545     $           (wm(i)*auxjn2(ind+1) + wp(i)*auxjn2(ind)) *
    549546     $           cortemp(i)           
    550547            !CO interpolated coefficient
    551548            jfotsout(indexint,11,auxi) =
    552      $           jfotsout(indexint,11,klev) *
     549     $           jfotsout(indexint,11,nlayer) *
    553550     $           (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind)) *
    554551     $           cortemp(i)           
    … …  
    556553         !Only if chemthermod.ge.2
    557554         if(chemthermod.ge.2) then
    558             do i=1,klev
     555            do i=1,nlayer
    559556               ind=auxind(i)
    560                auxi = klev-i+1
     557               auxi = nlayer-i+1
    561558               !NO interpolated coefficient
    562559               jfotsout(indexint,10,auxi)=
    563      $              jfotsout(indexint,10,klev) *
     560     $              jfotsout(indexint,10,nlayer) *
    564561     $              (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) *
    565562     $              cortemp(i)
    566563               !NO2 interpolated coefficient
    567564               jfotsout(indexint,13,auxi)=
    568      $              jfotsout(indexint,13,klev) *
     565     $              jfotsout(indexint,13,nlayer) *
    569566     $              (wm(i)*auxjno2(ind+1)+ wp(i)*auxjno2(ind)) *
    570567     $              cortemp(i)
    … …  
    585582c     Input atmospheric column
    586583
    587       do i=1,klev
    588          auxcolinp(klev-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
     584      do i=1,nlayer
     585         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
    589586     $        h2o2colx(i) + nocolx(i) + cocolx(i) + no2colx(i)
    590587      end do
    … …  
    620617         endif
    621618         call interfast
    622      $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
    623          do i=1,klev
     619     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
     620         do i=1,nlayer
    624621            ind=auxind(i)
    625             auxi = klev-i+1
     622            auxi = nlayer-i+1
    626623            !Correction to include T variation of CO2 cross section
    627624            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
    … …  
    634631            !CO2 interpolated coefficient
    635632            jfotsout(indexint,1,auxi) =
    636      $           jfotsout(indexint,1,klev) *
     633     $           jfotsout(indexint,1,nlayer) *
    637634     $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) *
    638635     $           cortemp(i) *
    … …  
    641638            !O2 interpolated coefficient
    642639            jfotsout(indexint,2,auxi) =
    643      $           jfotsout(indexint,2,klev) *
     640     $           jfotsout(indexint,2,nlayer) *
    644641     $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) *
    645642     $           cortemp(i)
    646643            !H2O interpolated coefficient
    647644            jfotsout(indexint,4,auxi) =
    648      $           jfotsout(indexint,4,klev) *
     645     $           jfotsout(indexint,4,nlayer) *
    649646     $           (wm(i)*auxjh2o(ind+1) + wp(i)*auxjh2o(ind)) *
    650647     $           cortemp(i)
    651648            !H2O2 interpolated coefficient
    652649            jfotsout(indexint,6,auxi) =
    653      $           jfotsout(indexint,6,klev) *
     650     $           jfotsout(indexint,6,nlayer) *
    654651     $           (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) *
    655652     $           cortemp(i)           
    656653            !CO interpolated coefficient
    657654            jfotsout(indexint,11,auxi) =
    658      $           jfotsout(indexint,11,klev) *
     655     $           jfotsout(indexint,11,nlayer) *
    659656     $           (wm(i)*auxjco(ind+1) + wp(i)*auxjco(ind)) *
    660657     $           cortemp(i)
    … …  
    662659         !Only if chemthermod.ge.2
    663660         if(chemthermod.ge.2) then
    664             do i=1,klev
     661            do i=1,nlayer
    665662               ind=auxind(i)
    666                auxi = klev-i+1
     663               auxi = nlayer-i+1
    667664               !NO interpolated coefficient
    668665               jfotsout(indexint,10,auxi)=
    669      $              jfotsout(indexint,10,klev) *
     666     $              jfotsout(indexint,10,nlayer) *
    670667     $              (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) *
    671668     $              cortemp(i)
    672669               !NO2 interpolated coefficient
    673670               jfotsout(indexint,13,auxi)=
    674      $              jfotsout(indexint,13,klev) *
     671     $              jfotsout(indexint,13,nlayer) *
    675672     $              (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind)) *
    676673     $              cortemp(i)
    … …  
    693690c     Input atmospheric column
    694691
    695       do i=1,klev
    696          auxcolinp(klev-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
     692      do i=1,nlayer
     693         auxcolinp(nlayer-i+1) = co2colx(i) + o2colx(i) + h2ocolx(i) +
    697694     $        h2o2colx(i) + nocolx(i) + no2colx(i)
    698695      end do
    … …  
    727724
    728725         call interfast
    729      $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
    730          do i=1,klev
     726     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
     727         do i=1,nlayer
    731728            ind=auxind(i)
    732             auxi = klev-i+1
     729            auxi = nlayer-i+1
    733730            !Correction to include T variation of CO2 cross section
    734731            if(sigma(indexint,auxi)*alfa(indexint,auxi)*
    … …  
    741738            !CO2 interpolated coefficient
    742739            jfotsout(indexint,1,auxi) =
    743      $           jfotsout(indexint,1,klev) *
     740     $           jfotsout(indexint,1,nlayer) *
    744741     $           (wm(i)*auxjco2(ind+1) + wp(i)*auxjco2(ind)) *
    745742     $           cortemp(i) *
    … …  
    748745            !O2 interpolated coefficient
    749746            jfotsout(indexint,2,auxi) =
    750      $           jfotsout(indexint,2,klev) *
     747     $           jfotsout(indexint,2,nlayer) *
    751748     $           (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) *
    752749     $           cortemp(i)
    753750            !H2O interpolated coefficient
    754751            jfotsout(indexint,4,auxi) =
    755      $           jfotsout(indexint,4,klev) *
     752     $           jfotsout(indexint,4,nlayer) *
    756753     $           (wm(i)*auxjh2o(ind+1) + wp(i)*auxjh2o(ind)) *
    757754     $           cortemp(i)
    758755            !H2O2 interpolated coefficient
    759756            jfotsout(indexint,6,auxi) =
    760      $           jfotsout(indexint,6,klev) *
     757     $           jfotsout(indexint,6,nlayer) *
    761758     $           (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) *
    762759     $           cortemp(i)           
    … …  
    764761         !Only if chemthermod.ge.2
    765762         if(chemthermod.ge.2) then
    766             do i=1,klev
     763            do i=1,nlayer
    767764               ind=auxind(i)
    768                auxi = klev-i+1
     765               auxi = nlayer-i+1
    769766               !NO interpolated coefficient
    770767               jfotsout(indexint,10,auxi)=
    771      $              jfotsout(indexint,10,klev) *
     768     $              jfotsout(indexint,10,nlayer) *
    772769     $              (wm(i)*auxjno(ind+1) +wp(i)*auxjno(ind)) *
    773770     $              cortemp(i)
    … …  
    795792
    796793      indexint=32
    797       do i=1,klev
    798          auxcolinp(klev-i+1) =co2colx(i) + o2colx(i) + h2o2colx(i) +
     794      do i=1,nlayer
     795         auxcolinp(nlayer-i+1) =co2colx(i) + o2colx(i) + h2o2colx(i) +
    799796     $        nocolx(i) + no2colx(i)
    800797      end do
    … …  
    824821      endif
    825822      call interfast
    826      $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
    827       do i=1,klev
     823     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
     824      do i=1,nlayer
    828825         ind=auxind(i)
    829          auxi = klev-i+1
     826         auxi = nlayer-i+1
    830827         !Correction to include T variation of CO2 cross section
    831          if(sigma(indexint,klev-i+1)*alfa(indexint,auxi)*
     828         if(sigma(indexint,nlayer-i+1)*alfa(indexint,auxi)*
    832829     $        coltemp(auxi).lt.60.) then
    833830            cortemp(i)=exp(-sigma(indexint,auxi)*
    … …  
    838835         !CO2 interpolated coefficient
    839836         jfotsout(indexint,1,auxi) =
    840      $        jfotsout(indexint,1,klev) *
     837     $        jfotsout(indexint,1,nlayer) *
    841838     $        (wm(i)*auxjco2(ind+1)+wp(i)*auxjco2(ind)) *
    842839     $        cortemp(i) *
    … …  
    845842         !O2 interpolated coefficient
    846843         jfotsout(indexint,2,auxi) =
    847      $        jfotsout(indexint,2,klev) *
     844     $        jfotsout(indexint,2,nlayer) *
    848845     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind)) *
    849846     $        cortemp(i)
    850847         !H2O2 interpolated coefficient
    851848         jfotsout(indexint,6,auxi) =
    852      $        jfotsout(indexint,6,klev) *
     849     $        jfotsout(indexint,6,nlayer) *
    853850     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind)) *
    854851     $        cortemp(i)         
    … …  
    856853      !Only if chemthermod.ge.2
    857854      if(chemthermod.ge.2) then
    858          do i=1,klev
    859             auxi = klev-i+1
     855         do i=1,nlayer
     856            auxi = nlayer-i+1
    860857            ind=auxind(i)
    861858            !NO interpolated coefficient
    862859            jfotsout(indexint,10,auxi) =
    863      $           jfotsout(indexint,10,klev) *
     860     $           jfotsout(indexint,10,nlayer) *
    864861     $           (wm(i)*auxjno(ind+1) + wp(i)*auxjno(ind)) *
    865862     $           cortemp(i)
    866863           !NO2 interpolated coefficient
    867864            jfotsout(indexint,13,auxi) =
    868      $           jfotsout(indexint,13,klev) *
     865     $           jfotsout(indexint,13,nlayer) *
    869866     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind)) *
    870867     $           cortemp(i)
    … …  
    885882
    886883      indexint=33
    887       do i=1,klev
    888          auxcolinp(klev-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
     884      do i=1,nlayer
     885         auxcolinp(nlayer-i+1) = o2colx(i) + h2o2colx(i) + no2colx(i)
    889886      end do
    890887
    … …  
    908905      endif
    909906      call interfast
    910      $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
    911       do i=1,klev
     907     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
     908      do i=1,nlayer
    912909         ind=auxind(i)
    913          auxi = klev-i+1
     910         auxi = nlayer-i+1
    914911         !O2 interpolated coefficient
    915          jfotsout(indexint,2,auxi) = jfotsout(indexint,2,klev) *
     912         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
    916913     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
    917914         !H2O2 interpolated coefficient
    918          jfotsout(indexint,6,auxi) = jfotsout(indexint,6,klev) *
     915         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
    919916     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
    920917      enddo
    921918      !Only if chemthermod.ge.2
    922919      if(chemthermod.ge.2) then
    923          do i=1,klev
     920         do i=1,nlayer
    924921            ind=auxind(i)
    925922            !NO2 interpolated coefficient
    926             jfotsout(indexint,13,klev-i+1) =
    927      $           jfotsout(indexint,13,klev) *
     923            jfotsout(indexint,13,nlayer-i+1) =
     924     $           jfotsout(indexint,13,nlayer) *
    928925     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
    929926         enddo
    … …  
    943940
    944941      indexint=34
    945       do i=1,klev
    946          auxcolinp(klev-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
     942      do i=1,nlayer
     943         auxcolinp(nlayer-i+1) = h2o2colx(i) + o2colx(i) + o3colx(i) +
    947944     $        no2colx(i)
    948945      end do
    … …  
    969966      endif
    970967      call interfast
    971      $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
    972       do i=1,klev
     968     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
     969      do i=1,nlayer
    973970         ind=auxind(i)
    974          auxi = klev-i+1
     971         auxi = nlayer-i+1
    975972         !O2 interpolated coefficient
    976          jfotsout(indexint,2,auxi) = jfotsout(indexint,2,klev) *
     973         jfotsout(indexint,2,auxi) = jfotsout(indexint,2,nlayer) *
    977974     $        (wm(i)*auxjo2(ind+1) + wp(i)*auxjo2(ind))
    978975         !H2O2 interpolated coefficient
    979          jfotsout(indexint,6,auxi) = jfotsout(indexint,6,klev) *
     976         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
    980977     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
    981978         !O3 interpolated coefficient
    982          jfotsout(indexint,7,auxi) = jfotsout(indexint,7,klev) *
     979         jfotsout(indexint,7,auxi) = jfotsout(indexint,7,nlayer) *
    983980     $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
    984981      enddo
    985982      !Only if chemthermod.ge.2
    986983      if(chemthermod.ge.2) then
    987          do i=1,klev
     984         do i=1,nlayer
    988985            ind=auxind(i)
    989986            !NO2 interpolated coefficient
    990             jfotsout(indexint,13,klev-i+1) =
    991      $           jfotsout(indexint,13,klev) *
     987            jfotsout(indexint,13,nlayer-i+1) =
     988     $           jfotsout(indexint,13,nlayer) *
    992989     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
    993990         enddo
    … …  
    10071004
    10081005      indexint=35
    1009       do i=1,klev
    1010          auxcolinp(klev-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
     1006      do i=1,nlayer
     1007         auxcolinp(nlayer-i+1) = h2o2colx(i) + o3colx(i) + no2colx(i)
    10111008      end do
    10121009
    … …  
    10301027      endif
    10311028      call interfast
    1032      $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
    1033       do i=1,klev
     1029     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
     1030      do i=1,nlayer
    10341031         ind=auxind(i)
    1035          auxi = klev-i+1
     1032         auxi = nlayer-i+1
    10361033         !H2O2 interpolated coefficient
    1037          jfotsout(indexint,6,auxi) = jfotsout(indexint,6,klev) *
     1034         jfotsout(indexint,6,auxi) = jfotsout(indexint,6,nlayer) *
    10381035     $        (wm(i)*auxjh2o2(ind+1) + wp(i)*auxjh2o2(ind))
    10391036         !O3 interpolated coefficient
    1040          jfotsout(indexint,7,auxi) = jfotsout(indexint,7,klev) *
     1037         jfotsout(indexint,7,auxi) = jfotsout(indexint,7,nlayer) *
    10411038     $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
    10421039      enddo
    10431040      if(chemthermod.ge.2) then
    1044          do i=1,klev
     1041         do i=1,nlayer
    10451042            ind=auxind(i)
    10461043            !NO2 interpolated coefficient
    1047             jfotsout(indexint,13,klev-i+1) =
    1048      $           jfotsout(indexint,13,klev) *
     1044            jfotsout(indexint,13,nlayer-i+1) =
     1045     $           jfotsout(indexint,13,nlayer) *
    10491046     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
    10501047         enddo
    … …  
    10631060
    10641061      indexint=36
    1065       do i=1,klev
    1066          auxcolinp(klev-i+1) = o3colx(i) + no2colx(i)
     1062      do i=1,nlayer
     1063         auxcolinp(nlayer-i+1) = o3colx(i) + no2colx(i)
    10671064      end do
    10681065
    … …  
    10841081      endif
    10851082      call interfast
    1086      $     (wm,wp,auxind,auxcolinp,klev,auxcoltab,nz2,limdown,limup)
    1087       do i=1,klev
     1083     $     (wm,wp,auxind,auxcolinp,nlayer,auxcoltab,nz2,limdown,limup)
     1084      do i=1,nlayer
    10881085         ind=auxind(i)
    10891086         !O3 interpolated coefficient
    1090          jfotsout(indexint,7,klev-i+1) =
    1091      $        jfotsout(indexint,7,klev) *
     1087         jfotsout(indexint,7,nlayer-i+1) =
     1088     $        jfotsout(indexint,7,nlayer) *
    10921089     $        (wm(i)*auxjo3(ind+1) + wp(i)*auxjo3(ind))
    10931090      enddo
    10941091      !Only if chemthermod.ge.2
    10951092      if(chemthermod.ge.2) then
    1096          do i=1,klev
     1093         do i=1,nlayer
    10971094            ind=auxind(i)
    10981095            !NO2 interpolated coefficient
    1099             jfotsout(indexint,13,klev-i+1) =
    1100      $           jfotsout(indexint,13,klev) *
     1096            jfotsout(indexint,13,nlayer-i+1) =
     1097     $           jfotsout(indexint,13,nlayer) *
    11011098     $           (wm(i)*auxjno2(ind+1) + wp(i)*auxjno2(ind))
    11021099         enddo
    … …  
    11071104c     End of interpolation to obtain photoabsorption rates
    11081105
    1109 
     1106c     Coefficients are refered to Sun-Mars distance
     1107c     Correction to the Sun-Venus distance (fixed)
     1108
     1109      jfotsout(:,:,:)=jfotsout(:,:,:)*(1.52/dist_sol)**2
     1110
     1111      end
     1112
     1113c**********************************************************************
     1114c**********************************************************************
     1115
     1116      subroutine column(ig,chemthermod,rm,nesptherm,tx,iz,zenit,
     1117     $     co2colx,o3pcolx, n2colx, cocolx)
     1118
     1119c     mar 2014        gg            adapted to Venus GCM
     1120c     nov 2002        fgg           first version
     1121
     1122c**********************************************************************
     1123      use dimphy
     1124      use param_v4_h
     1125      implicit none
     1126
     1127c    common variables and constants
     1128#include "clesphys.h"
     1129#include "mmol.h"
     1130
     1131c    local parameters and variables
     1132
     1133c     input and output variables
     1134
     1135      integer    ig
     1136      integer    chemthermod
     1137      integer    nesptherm                   !# of species undergoing
     1138chemistry, input
     1139      real       rm(klev,nesptherm)         !densities (cm-3), input
     1140      real       tx(klev)                   !temperature profile, input
     1141      real       iz(klev+1)                 !height profile, input
     1142      real       zenit                      !SZA, input
     1143      real       co2colx(klev)              !column density of CO2 (cm^-2), output
     1144      real       o3pcolx(klev)              !column density of O(3P)(cm^-2), output
     1145      real       n2colx(klev)               !N2 column density (cm-2), output
     1146      real       cocolx(klev)               !CO column density (cm-2), output
     1147
     1148c     local variables
     1149
     1150      real       xx
     1151      real       grav(klev)
     1152      real       Hco2,Ho3p,Ho2,Hh2,Hh2o,Hh2o2
     1153      real       Ho3,Hn2,Hn,Hno,Hco,Hh,Hno2
     1154
     1155      real       co2x(klev)
     1156      real       o3px(klev)
     1157      real       cox(klev)
     1158      real       n2x(klev)
     1159      real       nx(klev)
     1160
     1161      integer    i,j,k,icol,indexint          !indexes
     1162
     1163      integer    nz3
     1164
     1165      integer    jj
     1166      real*8      esp(klev*2)
     1167      real*8      ilayesp(klev*2)
     1168      real*8      szalayesp(klev*2)
     1169      integer     nlayesp
     1170      real*8      zmini
     1171      real*8      depth
     1172      real*8      espco2, espo2, espo3p, esph2, esph2o, esph2o2,espo3
     1173      real*8      espn2,espn,espno,espco,esph,espno2
     1174      real*8      rcmnz, rcmmini
     1175      real*8      szadeg
     1176
     1177! Tracer indexes in the thermospheric chemistry:
     1178!!! ATTENTION. These values have to be identical to those in euvheat.F90
     1179!!! If the values are changed there, the same has to be done here  !!!
     1180
     1181      integer,parameter :: ix_co2=1
     1182      integer,parameter :: ix_n2=13
     1183      integer,parameter :: ix_o=3
     1184      integer,parameter :: ix_co=4
     1185
     1186c*************************PROGRAM STARTS*******************************
     1187
     1188      nz3 = klev*2
     1189      do i=1,klev
     1190        xx = ( radio + iz(i) ) * 1.e5    ! conversion [km] ---> [cm]
     1191         grav(i) = gg * masa /(xx**2)    ! [cm/s2]
     1192      end do
     1193
     1194      !Scale heights  H = kT /Mg  -->  [cm]
     1195      xx = kboltzman * tx(klev) * n_avog / grav(klev)   ! g cm mol-1
     1196
     1197      Ho3p  = xx / mmolo
     1198      Hco2  = xx / mmolco2
     1199      Hco   = xx / mmolco
     1200      Hn2   = xx / mmoln2
     1201      Hn    = xx / mmoln
     1202
     1203      !Only if O3 chem. required
     1204c      if(chemthermod.ge.1)
     1205!     $     Ho3   = xx / mmol(igcm_o3)
     1206c     $     Ho3   = xx / mmolo3
     1207c      !Only if N or ion chem.
     1208c      if(chemthermod.ge.2) then
     1209c         Hn2   = xx / mmoln2
     1210c         Hn    = xx / mmoln
     1211c         Hno   = xx / mmolno
     1212c         Hno2  = xx / mmolno2
     1213c      endif
     1214      ! first loop in altitude : initialisation
     1215      do i=klev,1,-1
     1216         !Column initialisation
     1217         co2colx(i)  = 0.
     1218         o3pcolx(i)  = 0.
     1219         n2colx(i)   = 0.
     1220         cocolx(i)   = 0.
     1221
     1222         !--Densities [cm-3]
     1223         co2x(i)  = rm(i,ix_co2)
     1224         o3px(i)  = rm(i,ix_o)
     1225         cox(i)   = rm(i,ix_co)
     1226         n2x(i)   = rm(i,ix_n2)
     1227
     1228c         write(*,*), '--jthermcalc--', co2x(i)
     1229
     1230         !Only if O3 chem. required
     1231c         if(chemthermod.ge.1)
     1232c     $        o3x(i)   = rm(i,i_o3)
     1233c         !Only if Nitrogen of ion chemistry requested
     1234c         if(chemthermod.ge.2) then
     1235c            n2x(i)   = rm(i,i_n2)
     1236c            nx(i)    = rm(i,i_n)
     1237c            nox(i)   = rm(i,i_no)
     1238c            no2x(i)  = rm(i,i_no2)
     1239c         endif
     1240      enddo     ! end first loop
     1241
     1242      ! second loop in altitude : column calculations
     1243      do i=klev,1,-1
     1244         !Routine to calculate the geometrical length of each layer
     1245         call espesor_optico_A(ig,i,zenit,iz(i),nz3,iz,esp,ilayesp,
     1246     $         szalayesp,nlayesp, zmini)
     1247         if(ilayesp(nlayesp).eq.-1) then
     1248            co2colx(i)=1.e25
     1249            o3pcolx(i)=1.e25
     1250            n2colx(i)=1.e25
     1251            cocolx(i)=1.e25
     1252
     1253c            o2colx(i)=1.e25
     1254c            o3pcolx(i)=1.e25
     1255c            h2colx(i)=1.e25
     1256c            h2ocolx(i)=1.e25
     1257c            h2o2colx(i)=1.e25
     1258c            o3colx(i)=1.e25
     1259c            ncolx(i)=1.e25
     1260c            nocolx(i)=1.e25
     1261
     1262c            cocolx(i)=1.e25
     1263c            hcolx(i)=1.e25
     1264c            no2colx(i)=1.e25
     1265         else
     1266            rcmnz = ( radio + iz(klev) ) * 1.e5    ! km --> cm
     1267            rcmmini = ( radio + zmini ) * 1.e5
     1268            !Column calculation taking into account the geometrical
     1269            !depth
     1270            !calculated before
     1271            do j=1,nlayesp
     1272               jj=ilayesp(j)
     1273               !Top layer
     1274               if(jj.eq.klev) then
     1275                  if(zenit.le.60.) then
     1276                     o3pcolx(i)=o3pcolx(i)+o3px(klev)*Ho3p*esp(j)
     1277     $                    *1.e-5
     1278                     co2colx(i)=co2colx(i)+co2x(klev)*Hco2*esp(j)
     1279     $                    *1.e-5
     1280                     cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j)
     1281     $                    *1.e-5
     1282                     n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)
     1283     $                    *1.e-5
     1284
     1285c                     h2o2colx(i)=h2o2colx(i)+
     1286c     $                    h2o2x(klev)*Hh2o2*esp(j)*1.e-5
     1287c                     o2colx(i)=o2colx(i)+o2x(klev)*Ho2*esp(j)
     1288c     $                    *1.e-5
     1289c                     h2colx(i)=h2colx(i)+h2x(klev)*Hh2*esp(j)
     1290c     $                    *1.e-5
     1291c                     h2ocolx(i)=h2ocolx(i)+h2ox(klev)*Hh2o*esp(j)
     1292c     $                    *1.e-5                     
     1293c                     cocolx(i)=cocolx(i)+cox(klev)*Hco*esp(j)
     1294c     $                    *1.e-5
     1295c                     hcolx(i)=hcolx(i)+hx(klev)*Hh*esp(j)
     1296c     $                    *1.e-5
     1297                     !Only if O3 chemistry required
     1298c                     if(chemthermod.ge.1) o3colx(i)=
     1299c     $                    o3colx(i)+o3x(klev)*Ho3*esp(j)
     1300c     $                    *1.e-5
     1301                     !Only if N or ion chemistry requested
     1302c                     if(chemthermod.ge.2) then
     1303c                        n2colx(i)=n2colx(i)+n2x(klev)*Hn2*esp(j)
     1304c     $                    *1.e-5
     1305
     1306c                     endif
     1307                  else if(zenit.gt.60.) then
     1308                     espco2 =sqrt((rcmnz+Hco2)**2 -rcmmini**2) - esp(j)
     1309                     espo3p = sqrt((rcmnz+Ho3p)**2 -rcmmini**2)- esp(j)
     1310                     espco  = sqrt((rcmnz+Hco)**2 -rcmmini**2) - esp(j)
     1311                     espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
     1312                     espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
     1313
     1314c                     espo2  = sqrt((rcmnz+Ho2)**2 -rcmmini**2) - esp(j)
     1315c                     esph2  = sqrt((rcmnz+Hh2)**2 -rcmmini**2) - esp(j)
     1316c                     esph2o = sqrt((rcmnz+Hh2o)**2 -rcmmini**2)- esp(j)
     1317c                     esph2o2= sqrt((rcmnz+Hh2o2)**2-rcmmini**2)- esp(j)
     1318c                     esph   = sqrt((rcmnz+Hh)**2 -rcmmini**2)  - esp(j)
     1319                     !Only if O3 chemistry required
     1320c                     if(chemthermod.ge.1)                     
     1321c     $                   espo3=sqrt((rcmnz+Ho3)**2-rcmmini**2)-esp(j)
     1322c                     !Only if N or ion chemistry requested
     1323c                     if(chemthermod.ge.2) then
     1324c                        espn2 =sqrt((rcmnz+Hn2)**2-rcmmini**2)-esp(j)
     1325c                        espn  =sqrt((rcmnz+Hn)**2-rcmmini**2)  - esp(j)
     1326c                        espno =sqrt((rcmnz+Hno)**2-rcmmini**2) - esp(j)
     1327c                        espno2=sqrt((rcmnz+Hno2)**2-rcmmini**2)- esp(j)
     1328c                     endif
     1329
     1330                     co2colx(i) = co2colx(i) + espco2*co2x(klev)
     1331                     o3pcolx(i) = o3pcolx(i) + espo3p*o3px(klev)
     1332                     cocolx(i)  = cocolx(i)  + espco*cox(klev)
     1333                     n2colx(i)  = n2colx(i)  + espn2*n2x(klev)
     1334
     1335c                     o2colx(i)  = o2colx(i)  + espo2*o2x(klev)
     1336c                     h2colx(i)  = h2colx(i)  + esph2*h2x(klev)
     1337c                     h2ocolx(i) = h2ocolx(i) + esph2o*h2ox(klev)
     1338c                     h2o2colx(i)= h2o2colx(i)+ esph2o2*h2o2x(klev)
     1339c                     cocolx(i)  = cocolx(i)  + espco*cox(klev)
     1340c                     hcolx(i)   = hcolx(i)   + esph*hx(klev)
     1341                     !Only if O3 chemistry required
     1342c                     if(chemthermod.ge.1)                     
     1343c     $                  o3colx(i) = o3colx(i)  + espo3*o3x(klev)
     1344c                     !Only if N or ion chemistry requested
     1345c                     if(chemthermod.ge.2) then
     1346c                        n2colx(i)  = n2colx(i)  + espn2*n2x(klev)
     1347c                        ncolx(i)   = ncolx(i)   + espn*nx(klev)
     1348c                        nocolx(i)  = nocolx(i)  + espno*nox(klev)
     1349c                        no2colx(i) = no2colx(i) + espno2*no2x(klev)
     1350c                     endif
     1351                  endif !Of if zenit.lt.60
     1352               !Other layers
     1353               else
     1354                  co2colx(i)  = co2colx(i) +
     1355     $                 esp(j) * (co2x(jj)+co2x(jj+1)) / 2.
     1356                  o3pcolx(i)  = o3pcolx(i) +
     1357     $                 esp(j) * (o3px(jj)+o3px(jj+1)) / 2.
     1358                  cocolx(i)   = cocolx(i) +
     1359     $                 esp(j) * (cox(jj)+cox(jj+1)) / 2.
     1360                  n2colx(i)   = n2colx(i) +
     1361     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
     1362
     1363c
     1364c                  o2colx(i)   = o2colx(i) +
     1365c     $                 esp(j) * (o2x(jj)+o2x(jj+1)) / 2.
     1366c                  h2colx(i)   = h2colx(i) +
     1367c     $                 esp(j) * (h2x(jj)+h2x(jj+1)) / 2.
     1368c                  h2ocolx(i)  = h2ocolx(i) +
     1369c     $                 esp(j) * (h2ox(jj)+h2ox(jj+1)) / 2.
     1370c                  h2o2colx(i) = h2o2colx(i) +
     1371c     $                 esp(j) * (h2o2x(jj)+h2o2x(jj+1)) / 2.
     1372c                  hcolx(i)    = hcolx(i) +
     1373c     $                 esp(j) * (hx(jj)+hx(jj+1)) / 2.
     1374                  !Only if O3 chemistry required
     1375c                  if(chemthermod.ge.1)
     1376c     $                 o3colx(i) = o3colx(i) +
     1377c     $                 esp(j) * (o3x(jj)+o3x(jj+1)) / 2.
     1378c                  !Only if N or ion chemistry requested
     1379c                  if(chemthermod.ge.2) then
     1380c                     n2colx(i)   = n2colx(i) +
     1381c     $                 esp(j) * (n2x(jj)+n2x(jj+1)) / 2.
     1382c                     ncolx(i)    = ncolx(i) +
     1383c     $                    esp(j) * (nx(jj)+nx(jj+1)) / 2.
     1384c                     nocolx(i)   = nocolx(i) +
     1385c     $                    esp(j) * (nox(jj)+nox(jj+1)) / 2.
     1386c                     no2colx(i)  = no2colx(i) +
     1387c     $                    esp(j) * (no2x(jj)+no2x(jj+1)) / 2.
     1388c                  endif
     1389
     1390               endif  !Of if jj.eq.klev
     1391            end do    !Of do j=1,nlayesp
     1392         endif        !Of ilayesp(nlayesp).eq.-1
     1393      enddo           !Of do i=klev,1,-1
    11101394      return
    11111395
    … …  
    11131397
    11141398
    1115 
    1116 
    1117 
    1118 
     1399c**********************************************************************
     1400c**********************************************************************
     1401
     1402      subroutine interfast(wm,wp,nm,p,nlayer,pin,nl,limdown,limup)
     1403C
     1404C subroutine to perform linear interpolation in pressure from 1D profile
     1405C escin(nl) sampled on pressure grid pin(nl) to profile
     1406C escout(nlayer) on pressure grid p(nlayer).
     1407C
     1408      real*8 wm(nlayer),wp(nlayer),p(nlayer)
     1409      integer nm(nlayer)
     1410      real*8 pin(nl)
     1411      real*8 limup,limdown
     1412      integer nl,nlayer,n1,n,np,nini
     1413      nini=1
     1414      do n1=1,nlayer
     1415         if(p(n1) .gt. limup .or. p(n1) .lt. limdown) then
     1416            wm(n1) = 0.d0
     1417            wp(n1) = 0.d0
     1418         else
     1419            do n = nini,nl-1
     1420               if (p(n1).ge.pin(n).and.p(n1).le.pin(n+1)) then
     1421                  nm(n1)=n
     1422                  np=n+1
     1423                  wm(n1)=abs(pin(n)-p(n1))/(pin(np)-pin(n))
     1424                  wp(n1)=1.d0 - wm(n1)
     1425                  nini = n
     1426                  exit
     1427               endif
     1428            enddo
     1429         endif
     1430      enddo
     1431      return
     1432      end
     1433
     1434c**********************************************************************
     1435c**********************************************************************
     1436
     1437      subroutine espesor_optico_A (ig,capa, szadeg,z,
     1438     @                   nz3,iz,esp,ilayesp,szalayesp,nlayesp, zmini)
     1439
     1440c     fgg              nov 03      Adaptation to Martian model
     1441c     malv             jul 03      Corrected z grid. Split in alt & frec
     1442codes
     1443c     fgg              feb 03      first version
     1444*************************************************************************
     1445      use dimphy
     1446      use param_v4_h
     1447      implicit none
     1448
     1449c     arguments
     1450
     1451      real        szadeg                ! I. SZA [rad]
     1452      real        z                     ! I. altitude of interest [km]
     1453      integer     nz3,ig                   ! I. dimension of esp, ylayesp, etc...
     1454                                        !  (=2*klev= max# of layers in
     1455                                        !  ray path)
     1456      real     iz(klev+1)              ! I. Altitude of each layer
     1457      real*8        esp(nz3)            ! O. layer widths after geometrically
     1458                                        !    amplified; in [cm] except
     1459                                        !    at TOA
     1460                                        !    where an auxiliary value is
     1461                                        !    used
     1462      real*8        ilayesp(nz3)        ! O. Indexes of layers along ray path
     1463      real*8        szalayesp(nz3)      ! O. SZA [deg]    "     "     "
     1464      integer       nlayesp
     1465!      real*8        nlayesp             ! O. # layers along ray path at
     1466!      this z
     1467      real*8        zmini               ! O. Minimum altitud of ray path [km]
     1468
     1469c     local variables and constants
     1470
     1471        integer     j,i,capa
     1472        integer     jmin                  ! index of min.altitude along ray path
     1473        real*8      szarad                ! SZA [deg]
     1474        real*8      zz
     1475        real*8      diz(klev+1)
     1476        real*8      rkmnz                 ! distance TOA to center of Planet [km]
     1477        real*8      rkmmini               ! distance zmini to center of P [km]
     1478        real*8      rkmj                  ! intermediate distance to C of P [km]
     1479
     1480c external function
     1481        external  grid_R8          ! Returns index of layer containing the altitude
     1482                                ! of interest, z; for example, if
     1483                                ! zkm(i)=z or zkm(i)<z<zkm(i+1) =>
     1484                                ! grid(z)=i
     1485        integer   grid_R8
     1486
     1487*************************************************************************     
     1488        szarad = dble(szadeg)*3.141592d0/180.d0
     1489        zz=dble(z)
     1490        do i=1,klev
     1491           diz(i)=dble(iz(i))
     1492        enddo
     1493        do j=1,nz3
     1494          esp(j) = 0.d0
     1495          szalayesp(j) = 777.d0
     1496          ilayesp(j) = 0
     1497        enddo
     1498        nlayesp = 0
     1499
     1500        ! First case: szadeg<60
     1501        ! The optical thickness will be given by  1/cos(sza)
     1502        ! We deal with 2 different regions:
     1503        !   1: First, all layers between z and ztop ("upper part of
     1504        !   ray")
     1505        !   2: Second, the layer at ztop
     1506        if(szadeg.lt.60.d0) then
     1507
     1508           zmini = zz
     1509           if(abs(zz-diz(klev)).lt.1.d-3) goto 1357
     1510           ! 1st Zone: Upper part of ray
     1511           !
     1512           do j=grid_R8(zz,diz,klev),klev-1
     1513             nlayesp = nlayesp + 1
     1514             ilayesp(nlayesp) = j
     1515             esp(nlayesp) = (diz(j+1)-diz(j)) / cos(szarad)        ! [km]
     1516             esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
     1517             szalayesp(nlayesp) = szadeg
     1518           end do
     1519
     1520           !
     1521           ! 2nd Zone: Top layer
     1522 1357      continue
     1523           nlayesp = nlayesp + 1
     1524           ilayesp(nlayesp) = klev
     1525           esp(nlayesp) = 1.d0 / cos(szarad)         ! aux. non-dimens. factor
     1526           szalayesp(nlayesp) = szadeg
     1527
     1528        ! Second case:  60 < szadeg < 90
     1529        ! The optical thickness is evaluated.
     1530        !    (the magnitude of the effect of not using cos(sza) is 3.e-5
     1531        !     for z=60km & sza=30, and 5e-4 for z=60km & sza=60,
     1532        !     approximately)
     1533        ! We deal with 2 different regions:
     1534        !   1: First, all layers between z and ztop ("upper part of
     1535        !   ray")
     1536        !   2: Second, the layer at ztop ("uppermost layer")
     1537        else if(szadeg.le.90.d0.and.szadeg.ge.60.d0) then
     1538
     1539           zmini=(radio+zz)*sin(szarad)-radio
     1540           rkmmini = radio + zmini
     1541
     1542           if(abs(zz-diz(klev)).lt.1.d-4) goto 1470
     1543
     1544           ! 1st Zone: Upper part of ray
     1545           !
     1546           do j=grid_R8(zz,diz,klev),klev-1
     1547              nlayesp = nlayesp + 1
     1548              ilayesp(nlayesp) = j
     1549              esp(nlayesp) =
     1550     #             sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
     1551     #             sqrt( (radio+diz(j))**2 - rkmmini**2 )           ! [km]
     1552              esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
     1553              rkmj = radio+diz(j)
     1554              szalayesp(nlayesp) = asin( rkmmini/rkmj )             ! [rad]
     1555              szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592
     1556! [deg]
     1557           end do
     1558 1470      continue
     1559           ! 2nd Zone:  Uppermost layer of ray.
     1560           !
     1561           nlayesp = nlayesp + 1
     1562           ilayesp(nlayesp) = klev
     1563           rkmnz = radio+diz(klev)
     1564           esp(nlayesp)  =  sqrt( rkmnz**2 - rkmmini**2 )       ! aux.factor[km]
     1565           esp(nlayesp)  =  esp(nlayesp) * 1.d5                 ! aux.f. [cm]
     1566           szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
     1567           szalayesp(nlayesp) = szalayesp(nlayesp) * 180.d0/3.141592! [deg]
     1568
     1569
     1570        ! Third case:  szadeg > 90
     1571        ! The optical thickness is evaluated.
     1572        ! We deal with 5 different regions:
     1573        !   1: all layers between z and ztop ("upper part of ray")
     1574        !   2: the layer at ztop ("uppermost layer")
     1575        !   3: the lowest layer, at zmini
     1576        !   4: the layers increasing from zmini to z (here SZA<90)
     1577        !   5: the layers decreasing from z to zmini (here SZA>90)
     1578        else if(szadeg.gt.90.d0) then
     1579
     1580           zmini=(radio+zz)*sin(szarad)-radio
     1581           rkmmini = radio + zmini
     1582
     1583           if(zmini.lt.diz(1)) then         ! Can see the sun?  No => esp(j)=inft
     1584             nlayesp = nlayesp + 1
     1585             ilayesp(nlayesp) = - 1     ! Value to mark "no sun on view"
     1586!             esp(nlayesp) = 1.e30
     1587
     1588           else
     1589              jmin=grid_R8(zmini,diz,klev)+1
     1590             
     1591
     1592              if(abs(zz-diz(klev)).lt.1.d-4) goto 9876
     1593
     1594              ! 1st Zone: Upper part of ray
     1595              !
     1596              do j=grid_R8(zz,diz,klev),klev-1
     1597                nlayesp = nlayesp + 1
     1598                ilayesp(nlayesp) = j
     1599                esp(nlayesp) =
     1600     $                sqrt( (radio+diz(j+1))**2 - rkmmini**2 ) -
     1601     $                sqrt( (radio+diz(j))**2 - rkmmini**2 )          ! [km]
     1602                esp(nlayesp) = esp(nlayesp) * 1.d5                    ! [cm]
     1603                rkmj = radio+diz(j)
     1604                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
     1605                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592
     1606! [deg]
     1607              end do
     1608
     1609 9876         continue
     1610              ! 2nd Zone:  Uppermost layer of ray.
     1611              !
     1612              nlayesp = nlayesp + 1
     1613              ilayesp(nlayesp) = klev
     1614              rkmnz = radio+diz(klev)
     1615              esp(nlayesp) =  sqrt( rkmnz**2 - rkmmini**2 ) !aux.factor[km]
     1616              esp(nlayesp) = esp(nlayesp) * 1.d5 !aux.f.[cm]
     1617              szalayesp(nlayesp) = asin( rkmmini/rkmnz )           ! [rad]
     1618              szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592 ! [deg]
     1619
     1620              ! 3er Zone: Lowestmost layer of ray
     1621              !
     1622              if ( jmin .ge. 2 ) then      ! If above the planet's surface
     1623                j=jmin-1
     1624                nlayesp = nlayesp + 1
     1625                ilayesp(nlayesp) = j
     1626                esp(nlayesp) = 2. *
     1627     $                 sqrt( (radio+diz(j+1))**2 -rkmmini**2 )       ! [km]
     1628                esp(nlayesp) = esp(nlayesp) * 1.d5                   ! [cm]
     1629                rkmj = radio+diz(j+1)
     1630                szalayesp(nlayesp) = asin( rkmmini/rkmj ) ! [rad]
     1631                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592
     1632! [deg]
     1633              endif
     1634
     1635              ! 4th zone: Lower part of ray, increasing from zmin to z
     1636              !    ( layers with SZA < 90 deg )
     1637              do j=jmin,grid_R8(zz,diz,klev)-1
     1638                nlayesp = nlayesp + 1
     1639                ilayesp(nlayesp) = j
     1640                esp(nlayesp) =
     1641     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
     1642     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )       ! [km]
     1643                esp(nlayesp) = esp(nlayesp) * 1.d5                     ! [cm]
     1644                rkmj = radio+diz(j)
     1645                szalayesp(nlayesp) = asin( rkmmini/rkmj )              ! [rad]
     1646                szalayesp(nlayesp) = szalayesp(nlayesp) *180.d0/3.141592
     1647! [deg]
     1648              end do
     1649
     1650              ! 5th zone: Lower part of ray, decreasing from z to zmin
     1651              !    ( layers with SZA > 90 deg )
     1652              do j=grid_R8(zz,diz,klev)-1, jmin, -1
     1653                nlayesp = nlayesp + 1
     1654                ilayesp(nlayesp) = j
     1655                esp(nlayesp) =
     1656     $                    sqrt( (radio+diz(j+1))**2 - rkmmini**2 )
     1657     $                  - sqrt( (radio+diz(j))**2 - rkmmini**2 )
     1658! [km]
     1659                esp(nlayesp) = esp(nlayesp) * 1.d5
     1660! [cm]
     1661                rkmj = radio+diz(j)
     1662                szalayesp(nlayesp) = 3.141592 - asin( rkmmini/rkmj )
     1663! [rad]
     1664                szalayesp(nlayesp) = szalayesp(nlayesp)*180.d0/3.141592
     1665! [deg]
     1666              end do
     1667
     1668           end if
     1669
     1670        end if
     1671
     1672        return
     1673
     1674        end
     1675
     1676
     1677c**********************************************************************
     1678c***********************************************************************
     1679
     1680        function grid_R8 (z, zgrid, nz)
     1681
     1682c Returns the index where z is located within vector zgrid
     1683c The vector zgrid must be monotonously increasing, otherwise program stops.
     1684c If z is outside zgrid limits, or zgrid dimension is nz<2, the program stops.
     1685c
     1686c FGG     Aug-2004     Correct z.lt.zgrid(i) to .le.
     1687c MALV    Jul-2003
     1688c***********************************************************************
     1689
     1690        implicit none
     1691
     1692c Arguments
     1693        integer   nz
     1694        real*8      z
     1695        real*8      zgrid(nz)
     1696        integer   grid_R8
     1697
     1698c Local 
     1699        integer   i, nz1, nznew
     1700
     1701c*** CODE START
     1702
     1703        if ( z .lt. zgrid(1) .or. z.gt.zgrid(nz) ) then
     1704           write (*,*) ' GRID/ z outside bounds of zgrid '
     1705           write (*,*) ' z,zgrid(1),zgrid(nz) =', z,zgrid(1),zgrid(nz)
     1706           stop ' Serious error in GRID.F '
     1707        endif
     1708        if ( nz .lt. 2 ) then
     1709           write (*,*) ' GRID/ zgrid needs 2 points at least ! '
     1710           stop ' Serious error in GRID.F '
     1711        endif
     1712        if ( zgrid(1) .ge. zgrid(nz) ) then
     1713           write (*,*) ' GRID/ zgrid must increase with index'
     1714           stop ' Serious error in GRID.F '
     1715        endif
     1716
     1717        nz1 = 1
     1718        nznew = nz/2
     1719        if ( z .gt. zgrid(nznew) ) then
     1720           nz1 = nznew
     1721           nznew = nz
     1722        endif
     1723        do i=nz1+1,nznew
     1724           if ( z. eq. zgrid(i) ) then
     1725              grid_R8=i
     1726              return
     1727              elseif ( z .le. zgrid(i) ) then
     1728              grid_R8 = i-1
     1729              return
     1730           endif
     1731        enddo
     1732        grid_R8 = nz
     1733        return
     1734
     1735        end
     1736

  • trunk/LMDZ.VENUS/libf/phyvenus/load_ksi.F

    r2036 r2464  
    4444      real   lambdamin(nnuve),lambdamax(nnuve) ! in microns
    4545      real   dlambda                  ! in microns
     46      logical upper
    4647
    4748      nlve = klev
    … …  
    8788        read(10,*)
    8889        read(10,*) m,Nb
    89 cc      GG modif below 
    90         if (nlve.le.78.and.m.ne.nlve) then
    91          write(*,*) 'This is subroutine load_ksi'
    92          print*,'Dimension problem between ksi.txt and nlve'
    93          print*,'N levels = ',m,nlve
    94          stop
     90CC vertical axis check
     91        upper=.false.
     92cc the number of levels read for the matrix (m) should be the number of levels of the GCM (nlve=klev)
     93        if (m.ne.nlve) then
     94CC When GCM axis has more levels than the matrix, we must be in the case
     95cc with upper atmosphere, ie matrix limited to 78 levels and upper levels set to zero
     96cc otherwise there is a problem
     97         if ((m.eq.78).and.(nlve.gt.78)) then
     98          upper=.true.
     99         else
     100          write(*,*) 'This is subroutine load_ksi'
     101          print*,'Dimension problem between ksi.txt and nlve'
     102          print*,'N levels = ',m,nlve
     103          stop
     104         endif
    95105        endif
     106cc wavelength check
    96107        if (Nb.ne.nnuve) then
    97108         write(*,*) 'This is subroutine load_ksi'
    … …  
    113124            read(10,format_lect) (ksive(i,j,band,mat),j=0,m+1) ! no unit
    114125         enddo                  ! i
     126cc case upper atmosphere: the level for space (m+1) should be raised to nlve+1
     127         if (upper) then
     128           do i=0,m
     129             ksive(i,klev+1,band,mat)=ksive(i,m+1,band,mat)
     130             ksive(i,m+1,band,mat)=0.
     131           enddo
     132             ksive(klev+1,:,band,mat)=ksive(m+1,:,band,mat)
     133             ksive(m+1,:,band,mat)=0.
     134         endif
    115135        enddo                     ! band
    116136c       print*,"Matrice ",mat," lue"

  • trunk/LMDZ.VENUS/libf/phyvenus/moldiff_red.F90

    r1675 r2464  
    11subroutine moldiff_red(ngrid,nlayer,nq,pplay,pplev,pt,pq,&
    2                        ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff)
     2                       ptimestep,pdteuv,pdtconduc,pdqdiff)
    33
    44USE chemparam_mod
    … …  
    88implicit none
    99
    10 !#include "dimphys.h"
    1110#include "comcstfi.h"
    12 !#include "callkeys.h"
    13 !#include "comdiurn.h"
    14 !#include "chimiedata.h"
    15 !#include "tracer.h"
    16 !#include "conc.h"
    1711#include "diffusion.h"
    1812
    … …  
    2620      real ptimestep
    2721      real pplay(ngrid,nlayer)
    28       real zzlay(ngrid,nlayer)
    2922      real pplev(ngrid,nlayer+1)
    3023      real pq(ngrid,nlayer,nq)
    … …  
    6760      real*8,dimension(:),allocatable,save ::  wi,Wad,Uthermal,Lambdaexo,Hspecie
    6861      real*8,dimension(:),allocatable,save :: Mtot1,Mtot2,Mraf1,Mraf2
    69       character(len=20),dimension(14) :: ListeDiff
     62      integer,parameter :: nb_esp_diff=15
     63      character(len=20),dimension(nb_esp_diff) :: ListeDiff
    7064!ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    7165!     tracer numbering in the molecular diffusion
    … …  
    130124        ListeDiff(13)='o3'
    131125        ListeDiff(14)='n'
     126        ListeDiff(15)='he'
    132127        i_esc_h=1000
    133128        i_esc_h2=1000
    … …  
    140135        do nn=1,nq
    141136
    142         do n=1,14
     137        do n=1,nb_esp_diff
    143138        if (ListeDiff(n) .eq. tname(nn)) then
    144139        indic_diff(nn)=1
    … …  
    624619        enddo
    625620        enddo
     621
     622! DEBUG
     623!       print*,"rapport MASSE: ",Mtot2(:)/Mtot1(:)
    626624
    627625! Check mass  conservation of each specie on column

  • trunk/LMDZ.VENUS/libf/phyvenus/nirco2abs.F

    r2198 r2464  
    189189                  if(pq(ig,l,ico2) .gt. 1.e-6) then
    190190                     oco2gcm=pq(ig,l,io)/pq(ig,l,ico2)
    191 
     191                     ! handle the rare cases when pq(ig,l,io)<0
     192                     if (pq(ig,l,io).lt.0) then
     193                       write(*,*) "nirco2abs: warning ig=",ig," l=",l,
     194     &                            " pq(ig,l,io)=",pq(ig,l,io)
     195                       oco2gcm=1.e6
     196                     endif
    192197                  else
    193198                     oco2gcm=1.e6

  • trunk/LMDZ.VENUS/libf/phyvenus/nlte_setup.F

    r1442 r2464  
    177177                                !! k19 & k20
    178178
    179       k20x = 3.d-12
     179c      k20x = 3.d-12
    180180c  TEST GG: double the values of Kvv as recently found by Sharma et al.2014
    181 c      k20x = 6.d-12   
     181      k20x = 6.d-12   
    182182c  TEST GG: use the minimum value of the experimental bracket's values [1-6]
    183183c      k20x = 1.d-12

  • trunk/LMDZ.VENUS/libf/phyvenus/param_read_e107.F

    r1530 r2464  
    11      subroutine param_read_e107
    22
    3         use dimphy
    4         use conc
     3      use dimphy
     4      use param_v4_h, only: jfotsout,crscabsi2,
     5     .    c1_16,c17_24,c25_29,c30_31,c32,c33,c34,c35,c36,
     6     .    co2crsc195,co2crsc295,t0,
     7     .    jabsifotsintpar,ninter,nz2,
     8     .    nabs,e107,date_e107,e107_tab,
     9     .    coefit0,coefit1,coefit2,coefit3,coefit4,
     10     .    efdisco2,efdiso2,efdish2o,
     11     .    efdish2o2,efdish2,efdiso3,
     12     .    efdiso,efdisn,efdish,
     13     .    efdisno,efdisn2,efdisno2,
     14     .    efdisco,efionco2,efiono2,efionn2,
     15     .    efionco,efiono3p,efionn,
     16     .    efionno,efionh,
     17     .    fluxtop,ct1,ct2,p1,p2
     18
    519      implicit none
    620
    721
    822c     common variables and constants
    9       include 'param.h'
    10       include 'param_v4.h'
    11 c      include 'datafile.h'
    1223      include "clesphys.h"
    1324 
    … …  
    1930      real       nada
    2031      character*13  filename
    21       character (len=100) :: datafile="/u/forget/WWW/datagcm/datafile"
     32      character (len=100) :: datafile="HIGHATM"
    2233     
    2334c*************************+PROGRAM STARTS**************************

  • trunk/LMDZ.VENUS/libf/phyvenus/photochemistry_venus.F90

    r2453 r2464  
    14571457   cicol  = 0.
    14581458
    1459    do i = 1,nztable - 1
     1459   do i = 1,nztable-1
    14601460      if (table_colair(i) < col(iz)) then
    14611461         cicol = (log(col(iz)) - log(table_colair(i)))           &

  • trunk/LMDZ.VENUS/libf/phyvenus/physiq_mod.F

    r2286 r2464  
    6969      USE chemparam_mod
    7070      USE conc
     71      USE param_v4_h
    7172      USE compo_hedin83_mod2
    7273!      use ieee_arithmetic
    … …  
    257258      EXTERNAL nlthermeq
    258259      EXTERNAL euvheat
    259       EXTERNAL param_read
    260260      EXTERNAL param_read_e107
    261261      EXTERNAL conduction
    … …  
    356356c
    357357      REAL :: tr_seri(klon,klev,nqmax)
     358      REAL :: tr_hedin(klon,klev,nqmax)
    358359      REAL :: d_tr(klon,klev,nqmax)
    359360
    … …  
    411412
    412413c======================================================================
     414c======================================================================
    413415c INITIALISATIONS
    414 c================
     416c======================================================================
    415417
    416418      modname = 'physiq'
    … …  
    437439      END IF
    438440     
     441c======================================================================
    439442c PREMIER APPEL SEULEMENT
    440443c========================
    … …  
    584587
    585588         if (callthermos) then
    586             if(solvarmod.eq.0) call param_read
    587             if(solvarmod.eq.1) call param_read_e107
     589            call fill_data_thermos
     590            call allocate_param_thermos(klev)
     591            call param_read_e107
    588592         endif
    589593
    … …  
    597601            rho(ig,j)=pplay(ig,j)*mmean(ig,j)*1e-3/(rnew(ig,j)*t(ig,j))
    598602          enddo
    599 c        stop
    600603
    601604        enddo 
    … …  
    749752
    750753      if(callthermos) then
    751          call concentrations2(pplay,t,d_t,qx,nqmax)
     754         call concentrations2(pplay,t,qx,nqmax)
    752755      endif
    753756
     757c========================
    754758      ENDIF ! debut
    755 
     759c========================
     760
     761c======================================================================
    756762! ------------------------------------------------------
    757763!   Initializations done at every physical timestep:
    … …  
    923929c====================================================================
    924930c Orbite et eclairement
    925 c====================================================================
     931c=======================
    926932
    927933c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
    … …  
    931937     
    932938      zlongi = 0.0
    933       dist   = 0.72  ! en UA
     939      dist   = 0.72333  ! en UA
    934940
    935941c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
    … …  
    952958      end if
    953959
     960c======================================================================
     961c FIN INITIALISATIONS
     962c======================================================================
     963c======================================================================
     964
     965c=======================================================
     966! CONDUCTION  and  MOLECULAR VISCOSITY
     967c=======================================================
     968
     969        d_t_conduc(:,:)=0.
     970        d_u_molvis(:,:)=0.
     971        d_v_molvis(:,:)=0.
     972
     973        IF (callthermos) THEN
     974
     975           tsurf(:)=t_seri(:,1)
     976           call conduction(klon, klev,pdtphys,
     977     $            pplay,paprs,t_seri,
     978     $            tsurf,zzlev,zzlay,d_t_conduc)
     979
     980            call molvis(klon, klev,pdtphys,
     981     $            pplay,paprs,t_seri,
     982     $            u,tsurf,zzlev,zzlay,d_u_molvis)
     983
     984            call molvis(klon, klev, pdtphys,
     985     $            pplay,paprs,t_seri,
     986     $            v,tsurf,zzlev,zzlay,d_v_molvis)
     987
     988            DO k=1,klev
     989               DO ig=1,klon
     990                   t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
     991                   u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
     992                   v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
     993               ENDDO
     994            ENDDO
     995        ENDIF
     996
     997c====================================================================
     998c    Compute mean mass, cp and R :
     999c------------------------------------
     1000
     1001      if(callthermos) then
     1002         call concentrations2(pplay,t_seri,tr_seri, nqmax)
     1003      endif
     1004
     1005
     1006c=======================================================
     1007! CHEMISTRY AND MICROPHYSICS
     1008c=======================================================
     1009
    9541010      if (iflag_trac.eq.1) then
    9551011!====================================================================
    9561012! Case 1: pseudo-chemistry with relaxation toward fixed profile
    957 !====================================================================
     1013!=========
    9581014       if (tr_scheme.eq.1) then
    9591015
    … …  
    9681024! However, the variable 'source' could be used in physiq
    9691025! so the call to phytrac_emiss could be to initialise it.
    970 !====================================================================
    971 
    972          call phytrac_emiss ( (rjourvrai+gmtime)*RDAY,
    973      I                   debut,lafin,nqmax,
     1026!=========
     1027
     1028         call phytrac_emiss (debut,lafin,nqmax,
    9741029     I                   nlon,nlev,dtime,paprs,
    9751030     I                   latitude_deg,longitude_deg,
    … …  
    9851040!         nbapp_chem = 24000 => chempas = 4 => zctime = 420 s
    9861041!         nbapp_chem = 12000 => chempas = 8 => zctime = 840 s
    987 !====================================================================
     1042!=========
    9881043
    9891044         nbapp_chem = 24000
    … …  
    11311186
    11321187            do iq = 1, nqmax - nmicro
    1133                tr_seri(:,:,iq) = tr_seri(:,:,iq)
    1134      $                         + d_tr_chem(:,:,iq)*zctime
     1188               tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
     1189     $                         + d_tr_chem(:,:,iq)*zctime,1.e-30)
    11351190            end do
    11361191
    … …  
    11441199            else if (cl_scheme == 2) then
    11451200               do iq = nqmax-nmicro+1,nqmax
    1146                   tr_seri(:,:,iq) = tr_seri(:,:,iq)
    1147      $                            + d_tr_sed(:,:,iq)*zctime
     1201                  tr_seri(:,:,iq) = max(tr_seri(:,:,iq)
     1202     $                            + d_tr_sed(:,:,iq)*zctime,1.e-30)
    11481203               end do
    11491204            end if  ! cl_scheme
    … …  
    11521207         end if     ! mod(itap,chempas)  <------- end of chemistry supercycling
    11531208
    1154 !====================================================================
     1209!=========
    11551210! End Case 3: Full chemistry and/or clouds.
    11561211!====================================================================
    … …  
    13001355!      call writefield_phy('physiq_d_ts',d_ts,1)
    13011356
    1302 c
    1303 c Appeler l'ajustement sec
    1304 c
    13051357c===================================================================
    13061358c Convection seche
    … …  
    13711423      END IF
    13721424
    1373 c------------------------------------
    1374 c    Compute mean mass, cp and R :
    1375 c------------------------------------
    1376 
    1377       if(callthermos) then
    1378          call concentrations2(pplay,t_seri,d_t,tr_seri, nqmax)
    1379       endif
    1380 
    13811425c====================================================================
    13821426c RAYONNEMENT
    … …  
    13851429
    13861430      dtimerad = dtime*REAL(radpas)  ! pas de temps du rayonnement (s)
    1387 c====================================================================
    1388 
     1431
     1432! update mmean
     1433         if (ok_chem) then
     1434            mmean(:,:) = 0.
     1435            do iq = 1,nqmax - nmicro
     1436               mmean(:,:) = mmean(:,:)+tr_seri(:,:,iq)*m_tr(iq)
     1437               rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
     1438            enddo
     1439         endif
     1440
     1441cc---------------------------------------------
    13891442      if (callnlte .or. callthermos) then
    13901443         if (ok_chem) then
    … …  
    14111464        IF(callnlte) THEN                                 
    14121465            if(nltemodel.eq.0.or.nltemodel.eq.1) then
    1413                 CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
    1414      $                    tr_seri, d_t_nlte)
     1466c nltecool call not correct...
     1467c                CALL nltecool(klon, klev, nqmax, pplay*9.869e-6, t_seri,
     1468c     $                    tr_seri, d_t_nlte)
     1469                abort_message='nltemodel=0 or 1 should not be used...'
     1470                call abort_physic(modname,abort_message,1)
    14151471            else if(nltemodel.eq.2) then                               
     1472!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1473!! HEDIN instead of compo for this calculation
     1474!              call compo_hedin83_mod(pplay,rmu0,   
     1475!    $                 co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)     
     1476!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    14161477               CALL nlte_tcool(klon,klev,pplay*9.869e-6,             
    14171478     $               t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm,
    … …  
    14361497     $        CALL nlthermeq(klon, klev, paprs, pplay)
    14371498
    1438 c
     1499cc---------------------------------------------
    14391500c       LTE radiative transfert / solar / IR matrix
    14401501c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    14551516c       heat(:,:)=heat(:,:)/(1.+0.80*((rjourvrai-356)+gmtime)/20.)
    14561517
     1518cc---------------------------------------------
    14571519c       CO2 near infrared absorption
    14581520c      ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    14601522        d_t_nirco2(:,:)=0.
    14611523        if (callnirco2) then
     1524!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1525!! HEDIN instead of compo for this calculation
     1526!          call compo_hedin83_mod(pplay,rmu0,   
     1527!    $                 co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
     1528!          tr_hedin(:,:,:)=tr_seri(:,:,:)
     1529!          tr_hedin(:,:,i_co2)=co2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_co2)
     1530!          tr_hedin(:,:,i_co) = covmr_gcm(:,:)/mmean(:,:)*m_tr(i_co)
     1531!          tr_hedin(:,:,i_o)  =  ovmr_gcm(:,:)/mmean(:,:)*m_tr(i_o)
     1532!          tr_hedin(:,:,i_n2) = n2vmr_gcm(:,:)/mmean(:,:)*m_tr(i_n2)
     1533!          call nirco2abs (klon, klev, pplay, dist, nqmax, tr_hedin,
     1534!    .                 rmu0, fract, d_t_nirco2)
     1535!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    14621536           call nirco2abs (klon, klev, pplay, dist, nqmax, tr_seri,
    14631537     .                 rmu0, fract, d_t_nirco2)
     1538!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    14641539        endif
    14651540
    14661541
     1542cc---------------------------------------------
    14671543c          Net atmospheric radiative heating rate (K.s-1)
    14681544c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    14841560
    14851561cc---------------------------------------------
    1486 
    14871562c          EUV heating rate (K.s-1)
    14881563c          ~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    14921567        IF (callthermos) THEN
    14931568
    1494 c           call euvheat(klon, klev,t_seri,paprs,pplay,zzlay,
    1495 c     $          rmu0,pdtphys,gmtime,rjourvrai, co2vmr_gcm, n2vmr_gcm,
    1496 c     $          covmr_gcm, ovmr_gcm,d_t_euv )
    14971569           call euvheat(klon, klev, nqmax, t_seri,paprs,pplay,zzlay,
    1498      $         rmu0,dtimerad,gmtime,rjourvrai,
     1570     $         rmu0,dtimerad,gmtime,
     1571!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     1572!! HEDIN instead of compo for this calculation
     1573!! cf nlte_tcool and nirco2abs above !!
     1574!    $         tr_hedin, d_tr, d_t_euv )
     1575!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    14991576     $         tr_seri, d_tr, d_t_euv )
     1577!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    15001578                 
    15011579           DO k=1,klev
    … …  
    15081586        ENDIF  ! callthermos
    15091587
    1510 c====================================================================
    15111588      ENDIF    ! radpas
    15121589c====================================================================
    … …  
    15231600
    15241601      itap   = itap + 1
    1525 
    1526 ! CONDUCTION  and  MOLECULAR VISCOSITY
    1527 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1528 
    1529         d_t_conduc(:,:)=0.
    1530         d_u_molvis(:,:)=0.
    1531         d_v_molvis(:,:)=0.
    1532 
    1533         IF (callthermos) THEN
    1534 
    1535            tsurf(:)=t_seri(:,1)
    1536            call conduction(klon, klev,pdtphys,
    1537      $            pplay,paprs,t_seri,
    1538      $            tsurf,zzlev,zzlay,d_t_conduc)
    1539 
    1540             call molvis(klon, klev,pdtphys,
    1541      $            pplay,paprs,t_seri,
    1542      $            u,tsurf,zzlev,zzlay,d_u_molvis)
    1543 
    1544             call molvis(klon, klev, pdtphys,
    1545      $            pplay,paprs,t_seri,
    1546      $            v,tsurf,zzlev,zzlay,d_v_molvis)
    1547 
    1548             DO k=1,klev
    1549                DO ig=1,klon
    1550                   t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
    1551                   u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
    1552                   v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
    1553                ENDDO
    1554             ENDDO
    1555         ENDIF
    1556 
    1557 
     1602c====================================================================
     1603
     1604
     1605c==============================
    15581606!  --  MOLECULAR DIFFUSION ---
     1607c==============================
    15591608
    15601609          d_q_moldif(:,:,:)=0
    … …  
    15641613             call moldiff_red(klon, klev, nqmax,
    15651614     &                   pplay,paprs,t_seri, tr_seri, pdtphys,
    1566      &                   zzlay,d_t_euv,d_t_conduc,d_q_moldif)
     1615     &                   d_t_euv,d_t_conduc,d_q_moldif)
    15671616
    15681617
    … …  
    15721621           DO k=1,klev
    15731622              DO ig=1,klon
    1574                 tr_seri(ig,k,iq)= tr_seri(ig,k,iq)+
    1575      &                           d_q_moldif(ig,k,iq)*dtime ! [Kg/kg]?
     1623                tr_seri(ig,k,iq)= max(tr_seri(ig,k,iq)+
     1624     &                           d_q_moldif(ig,k,iq)*dtime,1.e-30)
    15761625              ENDDO
    15771626            ENDDO
    15781627           ENDDO
    15791628           
    1580 
    15811629         ENDIF  ! callthermos & ok_chem
    15821630
    … …  
    18781926      ENDDO
    18791927     
     1928c mise à jour rho,mmean pour sorties
     1929      if(callthermos) then
     1930         call concentrations2(pplay,t_seri,tr_seri, nqmax)
     1931      endif
     1932
    18801933c------------------------
    18811934c Calcul moment cinetique
    … …  
    20052058         do iq = 1,nqmax - nmicro
    20062059            call send_xios_field(tname(iq),
    2007      $                           qx(:,:,iq)*mmean(:,:)/m_tr(iq))
     2060     $                           tr_seri(:,:,iq)*mmean(:,:)/m_tr(iq))
    20082061         end do
    20092062
    … …  
    20122065         if ((tr_scheme == 3) .and. (cl_scheme == 1)) THEN
    20132066            call send_xios_field(tname(i_h2oliq),
    2014      $             qx(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq))
     2067     $             tr_seri(:,:,i_h2oliq)*mmean(:,:)/m_tr(i_h2oliq))
    20152068            call send_xios_field(tname(i_h2so4liq),
    2016      $             qx(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq))
     2069     $             tr_seri(:,:,i_h2so4liq)*mmean(:,:)/m_tr(i_h2so4liq))
    20172070            if (ok_sedim) then
    20182071               call send_xios_field("Fsedim",fsedim(:,1:klev))
    … …  
    20222075! chemical iterations
    20232076
    2024          call send_xios_field("iter",real(iter))
     2077         if (tr_scheme.eq.3) call send_xios_field("iter",real(iter))
    20252078
    20262079      end if
    … …  
    20312084       CALL send_xios_field("d_qmoldifN2",d_q_moldif(:,:,i_n2))
    20322085      ENDIF
     2086
     2087!! DEBUG
     2088!      if (is_master) print*,"itauphy=",itap
     2089!      if (itap.eq.10) lafin=.true.
    20332090
    20342091      if (lafin.and.is_omp_master) then

  • trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F

    r2320 r2464  
    1515
    1616      use chemparam_mod
    17       use conc, only: mmean
     17      use conc, only: mmean,rnew
    1818
    1919      implicit none
    … …  
    7373         if (reinit_trac .and. ok_chem) then
    7474
     75!!! in this reinitialisation, trac is VOLUME mixing ratio
    7576            print*, "Tracers are re-initialised"
    7677            trac(:,:,:) = 1.0e-30
    … …  
    106107            end if
    107108       
     109
     110! update mmean
     111            mmean(:,:) = 0.
     112            do iq = 1,nqmax - nmicro
     113               mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
     114               rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
     115            enddo
     116
    108117!           convert volume to mass mixing ratio
    109 
    110118            do iq = 1,nqmax - nmicro
    111119               trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:)
    … …  
    269277!===================================================================
    270278
     279! update mmean
     280            mmean(:,:) = 0.
     281            do iq = 1,nqmax - nmicro
     282               mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
     283               rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
     284            enddo
     285
    271286!     gas phase
    272287

  • trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F

    r2135 r2464  
    44c
    55c
    6       SUBROUTINE phytrac_emiss (timesimu,
    7      I                    debutphy,
     6      SUBROUTINE phytrac_emiss (debutphy,
    87     I                    lafin,
    98     I                    nqmax,
    109     I                    nlon,
    11      I                    nlev, 
     10     I                    nlev,
    1211     I                    pdtphys,
    13      I                    paprs,
     12     I                    paprs, 
    1413     I                    xlat,xlon,
    1514     O                    tr_seri)
     15     
     16     
    1617
    1718c======================================================================
    … …  
    2829c
    2930c======================================================================
    30       USE infotrac_phy, ONLY: nqtot
     31      use infotrac_phy, only: tname, nqtot
     32#ifdef CPP_XIOS     
     33      use xios_output_mod, only:  send_xios_field
     34#endif
    3135      use dimphy
    3236      USE geometry_mod, only: cell_area
    … …  
    3640#include "YOMCST.h"
    3741#include "clesphys.h"
     42     
    3843c======================================================================
    3944
    … …  
    4348c   ==========
    4449
    45       real timesimu   ! duree depuis debut simu (s)
    4650      logical debutphy       ! le flag de l'initialisation de la physique
    4751      logical lafin          ! le flag de la fin de la physique
    … …  
    6367cAA ----------------------------
    6468
    65 c pour emission volcan
     69c pour emission type effusion
    6670      real :: deltatr(klon,klev,nqtot)
    6771
    68       integer,parameter :: nblat=5,nblon=4,nbz=3
    69       integer,parameter :: Nemiss=0     ! duree emission (Ed)
    70       integer,save :: Nheight(nbz)      ! layer emission
    71       real,save :: so2_quantity         ! quantity so2 (kg)
    72       real,save :: lat_volcan(nblat),lon_volcan(nblon)
    73       real,save :: area_emiss(nblat,nblon)
    74       integer,save :: ig_volcan(nblat,nblon)
     72
     73
     74
     75      integer,parameter :: nblat=3,nblon=3,nbaire=3,nbflux=2,maxcell=16
     76      integer,parameter :: Nheight=3 ! layer emission (150m)
     77      integer,save :: Ncell(nbaire)
     78      real,save :: flux_surface_co2(nbflux) ! flux de CO2 emis (kg/m2/s)
     79      real :: flux(nlon,nqtot)
     80      real,save :: lat_zone(nblat),lon_zone(nblon)
     81      integer,save :: ig_zone(nblat,nblon,nbaire,nbflux,maxcell)
     82      integer,save :: numcell(nblat,nblon,nbaire,nbflux)
     83       
    7584
    7685      INTEGER i, k, it
    77       integer ilat,ilon,iz
     86      integer ilat,ilon,iaire,iflux,ipos
    7887      real    deltalat,deltalon
    79 c======================================================================
    80 
    81 c EMISSION TRACEURS
     88     
     89
     90
     91c======================================================================
     92
     93c EMISSIONS TRACEURS
    8294
    8395c---------
    … …  
    90102
    91103        ALLOCATE(M_tr(nqtot))
    92         M_tr(:)=64.                 ! SO2
    93        
     104        M_tr(:)=44.                ! CO2
     105
    94106C=========================================================================
    95107c Caracteristiques des traceurs emis:
    … …  
    97109
    98110c nombre total de traceur
    99          if (nbz*nblat*nblon .gt. nqtot) then
    100             print*, nbz*nblat*nblon, nqtot
     111         if (nblat*nblon*nbaire*nbflux+1 .gt. nqtot) then
     112            print*, nblat*nblon*nbaire*nbflux+1, nqtot
    101113            write(*,*) "Attention, pas assez de traceurs"
    102114            write(*,*) "le dernier sera bien le dernier"
    103          endif
    104 
    105 c quantite en kg
    106          so2_quantity = 20.*10.**9.
    107 
    108 c height (in layer index)
    109          Nheight(1) =  6  ! ~ 1 km
    110          Nheight(2) = 16  ! ~ 25 km
    111          Nheight(3) = 24  ! ~ 50 km
    112 
    113 c localisation volcan
    114          lat_volcan(1) =  70.
    115          lat_volcan(2) =  35.
    116          lat_volcan(3) =   0.
    117          lat_volcan(4) = -35.
    118          lat_volcan(5) = -70.
    119          lon_volcan(1) = -125.
    120          lon_volcan(2) =  -35.
    121          lon_volcan(3) =   55.
    122          lon_volcan(4) =  145.
    123          
    124          ig_volcan(ilat,ilon)= 0
     115         endif
     116                 
     117
     118
     119
     120c flux de CO2 (kg/s/m2)
     121         flux_surface_co2(1) = 5.*10.**-9.
     122         flux_surface_co2(2) = 5.*10.**-15.
     123
     124c nombre de cellule pour le cote du carre d'aire
     125         Ncell(1)= 2
     126         Ncell(2)= 3
     127         Ncell(3)= 4
     128     
     129
     130
     131c localisation zone emission
     132         lat_zone(1) = 08.
     133         lat_zone(2) = -50.
     134         lat_zone(3) = 35.
     135         lon_zone(1) = -172.
     136         lon_zone(2) = -20.
     137         lon_zone(3) = 70.
     138
     139 
    125140         if ((nbp_lon*nbp_lat)==1) then ! running a 1D simulation
    126141           deltalat=180.
    … …  
    131146         endif
    132147
     148         numcell(:,:,:,:)=0
     149         ig_zone(:,:,:,:,:)=0
    133150         do i=1,nlon
    134151          do ilat=1,nblat
    135152           do ilon=1,nblon
    136             if ((xlat(i).ge.lat_volcan(ilat))
    137      &     .and.((xlat(i)-deltalat).lt.lat_volcan(ilat))
    138      &     .and.(xlon(i).le.lon_volcan(ilon))
    139      &     .and.((xlon(i)+deltalon).gt.lon_volcan(ilon)) ) then
    140              ig_volcan(ilat,ilon)= i
    141              area_emiss(ilat,ilon) = cell_area(i)
    142              print*,"Lat,lon=",ilat,ilon," OK"
    143             end if
     153            do iaire=1,nbaire
     154
     155             if ((xlat(i).ge.lat_zone(ilat))
     156     &      .and.((xlat(i)-Ncell(iaire)*deltalat)
     157     &      .lt.lat_zone(ilat))
     158     &      .and.(xlon(i).le.lon_zone(ilon))
     159     &      .and.((xlon(i)+Ncell(iaire)*deltalon)
     160     &      .gt.lon_zone(ilon))) then
     161
     162              do iflux=1,nbflux
     163         numcell(ilat,ilon,iaire,iflux)=numcell(ilat,ilon,iaire,iflux)+1
     164        ig_zone(ilat,ilon,iaire,iflux,numcell(ilat,ilon,iaire,iflux))= i
     165             print*,"Lat,lon,naire,nflux,nlon=",ilat,ilon,iaire,iflux
     166     &        ,i," OK"
     167              end do
     168
     169             end if
     170
     171            end do
    144172           end do
    145173          end do
    … …  
    148176c Reinit des traceurs si necessaire
    149177         if (reinit_trac) then
    150            tr_seri(:,:,:)=0.
    151 c CAS N2 TRACEUR PASSIF POUR EVALUER PROFIL SOUS 7 KM
    152 c J'ai mis Nemiss=0 !
     178          tr_seri(:,:,:)=0.
     179
     180
    153181           do i=1,klon
    154182            do k=1,klev
    155               tr_seri(i,k,:)=max(min(0.035,
    156      &          0.035*(1.-log(paprs(i,k)/6.e6)/log(9.e6/6.e6))),0.)
    157             enddo
    158            enddo
    159 c FIN CAS N2 PASSIF
    160          endif
    161          
     183              tr_seri(i,k,:)=1.-28/43.44*max(min(0.035,
     184     &           0.035*(1.-log(paprs(i,k)/6.e6)/log(9.e6/6.e6))),0.)
     185            end do
     186           end do
     187         endif
     188 
    162189C=========================================================================
    163190C=========================================================================
    … …  
    168195
    169196c======================================================================
    170 c Emission d'un traceur pendant un certain temps
     197c Emission continue d'un traceur
    171198c necessite raz_date=1 dans run.def
    172199c et reinit_trac=y
    173200c======================================================================
    174201       deltatr(:,:,:) = 0.
    175 
    176 c source appliquee pendant Nemiss Ed
    177        if (timesimu .lt. 86400*Nemiss) then
     202       flux(:,:)=0.
    178203
    179204c emet les traceurs qui sont presents sur la grille
    180         do ilat  = 1,nblat
    181         do ilon  = 1,nblon
    182          if (ig_volcan(ilat,ilon).ne.0) then
    183          
    184           do iz = 1,nbz
    185            it=min( (iz-1)*nblat*nblon+(ilat-1)*nblon+ilon , nqtot )         
    186            i=ig_volcan(ilat,ilon)
    187            
     205      do ilat  = 1,nblat
     206       do ilon  = 1,nblon
     207        do iaire = 1,nbaire
     208         do iflux = 1,nbflux
     209     
     210              it=min( (ilat-1)*nblon*nbflux*nbaire+(iaire-1)*nbflux
     211     &         +(ilon-1)*nbaire*nbflux+iflux , nqtot )   
     212 
     213
    188214c injection dans une seule cellule:
    189215c source en kg/kg/s
    … …  
    196222
    197223c injection dans toute la colonne (a faire):
    198             do k=1,Nheight(iz)
    199              deltatr(i,k,it) = so2_quantity/(86400.*Nemiss) ! kg/s
    200      $  *RG/( area_emiss(ilat,ilon)
    201      $       *(paprs(i,1)-paprs(i,Nheight(iz)+1)) )    ! /kg (masse colonne)
    202      
    203              tr_seri(i,k,it) = tr_seri(i,k,it)+deltatr(i,k,it)*pdtphys
    204             end do
    205            
     224          do ipos=1,maxcell
     225 
     226           if (ig_zone(ilat,ilon,iaire,iflux,ipos).ne.0) then
     227
     228             i=ig_zone(ilat,ilon,iaire,iflux,ipos)         
     229             flux(i,it)=flux_surface_co2(iflux) 
     230             
     231            do k=1,Nheight
     232             deltatr(i,k,it) = flux_surface_co2(iflux) ! kg/s/m2
     233     $         *RG/(paprs(i,1)-paprs(i,Nheight+1))    ! /kg (masse colonne)
     234     
     235               tr_seri(i,k,it) = tr_seri(i,k,it)+deltatr(i,k,it)*pdtphys
     236            end do
     237
     238           end if
    206239          end do
    207          
    208          endif  ! ig_volcan!=0
    209         end do
    210         end do
    211 
    212        end if  ! duree emission
     240
     241        end do
     242        end do
     243        end do
     244        end do
     245
    213246       
    214247c======================================================================
    215248c======================================================================
     249
     250   
     251     
     252
     253     
     254#ifdef CPP_XIOS     
     255       do it=1,nqtot
     256       CALL  send_xios_field("flux_"//tname(it),
     257     &                     flux(:,it))
     258      end do
     259#endif
     260
     261
    216262
    217263      RETURN

  • trunk/LMDZ.VENUS/libf/phyvenus/sw_venus_rh.F

    r2048 r2464  
    285285c         PHEAT(j) = PHEAT(j)*2.
    286286c        endif
    287         if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
     287c BASE:
     288c        if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
     289c          PHEAT(j) = PHEAT(j)*3
     290c        endif
     291c AFTER Tayloring TdeepD
     292         if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
     293           PHEAT(j) = PHEAT(j)*3.5
     294         endif
     295         if ((PPB(j).gt.10.).and.(PPB(j).le.50.)) then
     296           PHEAT(j) = PHEAT(j)*1.5
     297         endif
     298c Options:
     299c        if ((PPB(j).gt.1.4).and.(PPB(j).le.10.)) then
     300c        if ((PPB(j).gt.2.0).and.(PPB(j).le.10.)) then
    288301c        if ((PPB(j).gt.1.4).and.(PPB(j).le.100.)) then
    289          PHEAT(j) = PHEAT(j)*3
    290         endif
    291 c        if ((PPB(j).gt.10.).and.(PPB(j).le.120.)) then
    292 c         PHEAT(j) = PHEAT(j)*2
     302c          PHEAT(j) = PHEAT(j)*3.5
     303c          PHEAT(j) = PHEAT(j)*3
     304c          PHEAT(j) = PHEAT(j)*2.5
     305c        endif
     306c        if ((PPB(j).gt.10.).and.(PPB(j).le.35.)) then
     307c        if ((PPB(j).gt.10.).and.(PPB(j).le.50.)) then
     308c          PHEAT(j) = PHEAT(j)*2
     309c          PHEAT(j) = PHEAT(j)*1.5
     310c          PHEAT(j) = PHEAT(j)*1.3
     311c        endif
     312c        if ((PPB(j).gt.35.).and.(PPB(j).le.120.)) then
     313c          PHEAT(j) = PHEAT(j)*2
     314c          PHEAT(j) = PHEAT(j)*1.5
    293315c        endif
    294316c----------------
Note: See TracChangeset for help on using the changeset viewer.