Ignore:
Timestamp:
Dec 17, 2013, 1:02:44 PM (11 years ago)
Author:
slebonnois
Message:

SL: update of Titan photochemical module to include computation of chemistry up to 1300 km

Location:
trunk/LMDZ.TITAN/libf/phytitan
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/phytitan/calchim.F

    r1056 r1126  
    11      SUBROUTINE calchim(nlon,ny,qy_c,nomqy_c,declin_rad,ls_rad,dtchim,
    2      .                   ctemp,cplay,cplev,
     2     .                   ctemp,cplay,cplev,czlay,czlev,
    33     .                   dqyc)
    44     
    … …  
    1010c            adaptation pour Titan 3D: 02/2009
    1111c            adaptation pour // : 04/2013
    12 c
     12c            extension chimie jusqu a 1300km : 10/2013
     13c
    1314c-------------------------------------------------
    1415c
    1516      use dimphy
    1617      use common_mod, only:utilaer,maer,prodaer,csn,csh,psurfhaze,
    17      .                     NLEV,NC,ND,NR
     18     .                     NLEV,NLRT,NC,ND,NR
    1819      USE comgeomphy,  only: rlatd
    19       use moyzon_mod,  only: klat
     20      USE moyzon_mod, only: tmoy,playmoy,zlaymoy,zlevmoy,klat
    2021      implicit none
    2122#include "dimensions.h"
    … …  
    3536      REAL         ctemp(nlon,klev)      ! Temperature
    3637      REAL         cplay(nlon,klev)      ! pression (Pa)
    37       REAL         cplev(nlon,klev)      ! pression intercouches (Pa)
     38      REAL         cplev(nlon,klev+1)    ! pression intercouches (Pa)
     39      REAL         czlay(nlon,klev)      ! altitude (m)
     40      REAL         czlev(nlon,klev+1)    ! altitude intercouches (m)
    3841     
    3942      REAL         dqyc(nlon,klev,NC)    ! Tendances especes chimiques
    … …  
    4649c variables envoyees dans la chimie: double precision
    4750
    48       REAL  temp_c(klev),press_c(klev)     ! T,p(mbar) a 1 lat donnee
    49       REAL  declin_c                       ! declinaison en degres
    50       REAL  cqy(klev,NC)                   ! frac mol qui seront modifiees
    51       REAL  surfhaze(klev)
    52       REAL  cprodaer(klev,4),cmaer(klev,4),ccsn(klev,4),ccsh(klev,4)
    53       REAL  rinter(klev),nb(klev)
     51      REAL  temp_c(NLEV)
     52      REAL  press_c(NLEV)   ! T,p(mbar) a 1 lat donnee
     53      REAL  cqy(NLEV,NC)    ! frac mol qui seront modifiees
     54      REAL  cqy0(NLEV,NC)    ! frac mol avant modif
     55      REAL  surfhaze(NLEV)
     56      REAL  cprodaer(NLEV,4),cmaer(NLEV,4)
     57      REAL  ccsn(NLEV,4),ccsh(NLEV,4)
     58! rmil: milieu de couche, grille pour K,D,p,ct,T,y
     59! rinter: intercouche (grille RA dans la chimie)
     60      REAL  rmil(NLEV),rinter(NLEV),nb(NLEV)
     61      REAL,save :: kedd(NLEV)
     62
    5463      REAL  a,b
    5564      character str1*1,str2*2
    56       integer ntotftop
    57       parameter (ntotftop=50)
    58       integer nftop(ntotftop),isaisonflux
    59       REAL  fluxtop(NC),latit,tabletmp(ntotftop),factflux
    60       character*10 nomsftop(ntotftop+1)
     65      REAL  latit
    6166      character*20 formt1,formt2,formt3
    6267     
    … …  
    6772      SAVE    firstcal
    6873     
    69       REAL  mass(NC),duree
    70       REAL  tablefluxtop(NC,jjp1,5)
    71       REAL  botCH4
     74      integer dec,declinint,ialt
     75      REAL  declin_c                       ! declinaison en degres
     76      real  factalt,factdec,krpddec,krpddecp1,krpddecm1
     77      real  duree
     78      REAL,save :: mass(NC)
     79      REAL,save,allocatable :: md(:,:)
     80      REAL,save :: botCH4
    7281      DATA  botCH4/0.05/
     82      real,save ::  r1d(131),ct1d(131),p1d(131),t1d(131) ! lecture tcp.ver
    7383      REAL,save,allocatable :: krpd(:,:,:,:),krate(:,:)
    74       integer reactif(5,NR),nom_prod(NC),nom_perte(NC)
    75       integer prod(200,NC),perte(2,200,NC)
    76       SAVE    mass,tablefluxtop,botCH4
    77       SAVE    reactif,nom_prod,nom_perte,prod,perte
    78      
     84      integer,save :: reactif(5,NR),nom_prod(NC),nom_perte(NC)
     85      integer,save :: prod(200,NC),perte(2,200,NC)
     86      character  dummy*30,fich*7,ficdec*15,curdec*15,name*10
     87      real  ficalt,oldq,newq,zalt
     88      logical okfic
     89
    7990c-----------------------------------------------------------------------
    8091c***********************************************************************
    … …  
    94105c ************************************
    95106
    96         allocate(krpd(15,ND+1,klev,jjp1),krate(klev,NR))
    97 
    98 c Verification dimension verticale: coherence titan_for.h et klev
    99 c --------------------------------
    100 
    101         if (klev.ne.NLEV) then
    102            print*,'PROBLEME de coherence dans la dimension verticale:'
    103      .           ,klev,NLEV
    104            STOP
    105         endif
    106 
    107 c Verification du nombre de composes: coherence titan_for.h et nqmax-nmicro
     107        allocate(krpd(15,ND+1,NLRT,jjp1),krate(NLEV,NR),md(NLEV,NC))
     108
     109c Verification du nombre de composes: coherence common_mod et nqmax-nmicro
    108110c ----------------------------------
    109111
    … …  
    114116        endif
    115117
    116 c calcul de temp_c, densites et press_c au milieu de l'ensemble des points:
    117 c ----------------------------------------------------------------------
    118 
    119         print*,'pression, densites et temp (chimie):'
     118c calcul de temp_c, densites et press_c en moyenne planetaire :
     119c --------------------------------------------------------------
     120
     121        print*,'pression, densites et temp (init chimie):'
    120122        print*,'level, press_c, nb, temp_c'
    121123        DO l=1,klev
     124         rinter(l)  = (zlevmoy(l)+RA)/1000.
     125         rmil(l)    = (zlaymoy(l)+RA)/1000.
    122126c     temp_c (K):
    123          temp_c(l)  = ctemp(nlon/2+1,l)
     127         temp_c(l)  = tmoy(l)
    124128c     press_c (mbar):
    125          press_c(l) = cplay(nlon/2+1,l)/100.
     129         press_c(l) = playmoy(l)/100.
    126130c     nb (cm-3):
    127131         nb(l) = 1.e-4*press_c(l) / (RKBOL*temp_c(l))
    128          print*,l,press_c(l),nb(l),temp_c(l)
     132         print*,l,rmil(l)-RA/1000.,press_c(l),nb(l),temp_c(l)
    129133        ENDDO
     134
     135c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
     136       do l=klev+1,NLEV
     137         rinter(l) = rinter(klev)
     138     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
     139         rmil(l)   = rmil(klev)
     140     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
     141       enddo       
     142
     143c lecture de tcp.ver, une seule fois
     144c remplissage r1d,t1d,ct1d,p1d
     145        open(11,file='../../INPUT/tcp.ver',status='old')
     146        read(11,*) dummy
     147        do i=1,131
     148          read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i)
     149c         print*,"TCP.VER ",r1d(i),t1d(i),ct1d(i),p1d(i)
     150        enddo
     151        close(11)
     152
     153c extension pour klev+1 a NLEV avec tcp.ver
     154c la jonction klev/klev+1 est brutale... Tant pis ?
     155        ialt=1
     156        do l=klev+1,NLEV
     157           zalt=rmil(l)-RA/1000.
     158           do i=ialt,130
     159            if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
     160              ialt=i
     161            endif
     162           enddo
     163           factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
     164           press_c(l) = exp(  log(p1d(ialt))   * (1-factalt)
     165     &                      + log(p1d(ialt+1)) * factalt    )
     166           nb(l)      = exp(  log(ct1d(ialt))   * (1-factalt)
     167     &                      + log(ct1d(ialt+1)) * factalt    )
     168           temp_c(l)  = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
     169           print*,l,zalt,press_c(l),nb(l),temp_c(l)
     170        enddo
    130171       
    131172c caracteristiques des composes:       
    132173c -----------------------------
    133         mass = 0.0
    134         call comp(nomqy_c,mass)
     174        mass(:) = 0.0
     175        call comp(nomqy_c,nb,temp_c,mass,md)
    135176        print*,'           Mass'
    136177        do ic=1,NC
    137178          print*,nomqy_c(ic),mass(ic)
     179c         if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,ic)
    138180        enddo
    139181       
    140182
    141 c FLUX AU SOMMET: DEPENDANT DE LA LATITUDE ET DE LA SAISON...
    142 
    143 c flux dans la derniere couche:(issu du modele 1D, a 500 km)
    144 c -----------------------------
    145 c       !!  en cm-2.s-1 !!
    146 c  ces valeurs sont converties en cm-3.s-1 (dni/dt) dans a couche 54
    147 c
    148 c Lecture des tables de flux en fonction lat et saison
    149 c ----------------------------------------------------
    150 
    151         WRITE(str2(1:2),'(i2.2)') ntotftop
    152         formt1 = '(A10,'//str2//'(3X,A10))'
    153         formt2 = '(F12.6,'//str2//'(2X,F13.6))'
    154         formt3 = '(F13.6,'//str2//'(2X,F13.6))'
    155        
    156         do j=1,jjp1
    157           do ic=1,NC
    158             do l=1,5
    159              tablefluxtop(ic,j,l) = 0.
    160             enddo
    161           enddo
    162         enddo
    163        
    164         do l=1,4
    165           WRITE(str1(1:1),'(i1.1)') l
    166 c hx1 -> Ls=RPI/2 / hx4 -> Ls=0
    167           open(10,file="../../INPUT/FLUXTOP/flux500.hs"//str1)
    168 c         open(10,file="FLUXTOP/flux500.hs"//str1)
    169 c on remet l=1 et 5 -> Ls=0 / l=2 -> Ls=RPI/2 ...
    170 
    171           read(10,formt1) nomsftop
    172           do j=1,ntotftop
    173            do i=1,10   !justification a gauche...
    174              if (nomsftop(j+1)(i:i).ne.' ') then
    175                 nomsftop(j) = nomsftop(j+1)(i:)
    176                 goto 100
    177              endif
    178            enddo
    179 100        continue
    180 c         print*,nomsftop(j)
    181            nftop(j) = 0
    182            do i=1,NC
    183             if (nomqy_c(i).eq.nomsftop(j)) then
    184              nftop(j) = i
    185             endif
    186            enddo
    187 c          if (nftop(j).ne.0) print*,nomsftop(j),nomqy_c(nftop(j))
    188           enddo
    189 
    190          
    191 c lecture des flux. Les formats sont un peu alambiques...
    192 c     a revoir eventuellement
    193           do j=1,jjp1/2+1
    194             read(10,formt2) latit,tabletmp
    195             do i=1,ntotftop
    196               if (nftop(i).ne.0) then
    197                tablefluxtop(nftop(i),j,l+1) = tabletmp(i)
    198               endif
    199             enddo
    200           enddo
    201           do j=jjp1/2+2,jjp1
    202             read(10,formt3) latit,tabletmp
    203             do i=1,ntotftop
    204               if (nftop(i).ne.0) then
    205                tablefluxtop(nftop(i),j,l+1) = tabletmp(i)
    206               endif
    207             enddo
    208           enddo
    209 
    210           close(10)
    211         enddo  ! l
    212        
    213         do j=1,jjp1
    214           do ic=1,NC
    215              tablefluxtop(ic,j,1) = tablefluxtop(ic,j,5)
    216           enddo
    217         enddo
    218        
    219 c  Stockage des composes utilises dans la prod d'aerosols
     183c  Stockage des composes utilises dans la prod d aerosols
    220184c     (aerprod=1) et dans H-> H2 (htoh2=1): utilaer
    221185c     ! decalage de 1 car utilise dans le c !
    222 c ------------------------------------------
    223 c + modifs du flux en fonction des obs UVIS et CIRS: C2H2,C2H6,HCN
    224 c + modifs des flux qui presentent un probleme evident: C3H8 et C4H10
    225 c ------------------------------------------
    226186
    227187        do ic=1,NC
    … …  
    244204          if (nomqy_c(ic).eq."HCN") then
    245205            utilaer(6)  = ic-1
    246             do j=1,jjp1
    247               do i=1,5
    248              tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 6.
    249               enddo
    250             enddo
    251206          endif
    252207          if (nomqy_c(ic).eq."HC3N") then
    … …  
    268223          if (nomqy_c(ic).eq."C2H2") then
    269224            utilaer(3)  = ic-1
    270             do j=1,jjp1
    271               do i=1,5
    272 c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 7.3
    273              tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 10.
    274               enddo
    275             enddo
    276225          endif
    277226          if (nomqy_c(ic).eq."AC6H6") then
    … …  
    298247          if (nomqy_c(ic).eq."AC6H5") then
    299248            utilaer(12) = ic-1
    300           endif
    301 
    302           if (nomqy_c(ic).eq."C2H6") then
    303             do j=1,jjp1
    304               do i=1,5
    305 c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 640.
    306              tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * 1200.
    307               enddo
    308             enddo
    309           endif
    310           if (nomqy_c(ic).eq."C3H8") then
    311             do j=1,jjp1
    312               do i=1,5
    313 c             tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.)
    314              tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-3000.)
    315               enddo
    316             enddo
    317           endif
    318           if (nomqy_c(ic).eq."C4H10") then
    319             do j=1,jjp1
    320               do i=1,5
    321              tablefluxtop(ic,j,i) = tablefluxtop(ic,j,i) * (-1.)
    322               enddo
    323             enddo
    324249          endif
    325250
    … …  
    369294c       enddo
    370295
     296
     297c   init kedd
     298c   ---------
     299c   kedd stays constant with time and space
     300c   => linked to pressure rather than altitude...
     301 
     302      kedd(:) = 5e2 ! valeur mise par defaut 
     303               ! UNITE ? doit etre ok pour gptitan
     304      do l=1,NLEV
     305       zalt=rmil(l)-RA/1000.  ! z en km
     306       if     ((zalt.ge.250.).and.(zalt.le.400.)) then
     307         kedd(l) = 10.**(3.+(zalt-250.)/50.)
     308! 1E3 at 250 km
     309! 1E6 at 400 km
     310       elseif ((zalt.gt.400.).and.(zalt.le.600.)) then
     311         kedd(l) = 10.**(6.+1.3*(zalt-400.)/200.)
     312! 2E7 at 600 km
     313       elseif ((zalt.gt.600.).and.(zalt.le.900.)) then
     314         kedd(l) = 10.**(7.3+0.7*(zalt-600.)/300.)
     315! 1E8 above 900 km
     316       elseif ( zalt.gt.900.                    ) then
     317        kedd(l) = 1.e8
     318       endif
     319      enddo
     320
    371321      ENDIF  ! premier appel
    372322
    … …  
    380330c      print*,'declinaison en degre=',declin_c
    381331       
    382 c-----------------------------------------------------------------------
    383 c
    384 c   calcul du facteur d'interpolation entre deux saisons pour flux
    385 c   --------------------------------------------------------------
    386 
    387        isaisonflux = int(ls_rad*2./RPI)+1
    388        if (isaisonflux.eq.5) then
    389          isaisonflux = 1
    390        endif
    391        factflux =    (sin(ls_rad)-sin((isaisonflux-1)*RPI/2.))/
    392      .   (sin(isaisonflux*RPI/2.)-sin((isaisonflux-1)*RPI/2.))
    393 
    394332c***********************************************************************
    395333c***********************************************************************
    … …  
    410348c***********************************************************************
    411349
     350c-----------------------------------------------------------------------
     351c
     352c   Distance radiale (intercouches, en km)
     353c   ----------------------------------------
     354
     355       do l=1,klev
     356         rinter(l) = (RA+czlev(j,l))/1000.
     357         rmil(l)   = (RA+czlay(j,l))/1000.
     358c        print*,rinter(l)
     359       enddo
     360
     361c au-dessus du GCM, dr regulier et rinter(NLEV)=1290+2575 km.
     362       do l=klev+1,NLEV
     363         rinter(l) = rinter(klev)
     364     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
     365         rmil(l)   = rmil(klev)
     366     &          + (l-klev)*(3865.-rinter(klev))/(NLEV-klev)
     367       enddo
     368       
    412369c-----------------------------------------------------------------------
    413370c
    … …  
    423380               nb(l) = 1.e-4*press_c(l) / (RKBOL*temp_c(l))
    424381       ENDDO
    425        
    426 c-----------------------------------------------------------------------
    427 c
    428 c   Distances radiales (intercouches, en km)
    429 c   ----------------------------------------
    430 
    431        rinter(1) = RA/1000. 
    432        do l=2,klev
    433 c A REVOIR: g doit varier avec r !
    434          rinter(l) = rinter(l-1) +
    435      .    (cplev(j,l-1)-cplev(j,l))/cplay(j,l)*(RD*temp_c(l)/RG)/1000.
    436 c        print*,rinter(l)
    437        enddo
    438 
    439 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    440 c FLUX AU TOP: interpolation en saison et
    441 c conversion en cm-3 s-1 dans la couche klev-1 (NLEV-2 dans le C)
    442 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    443         do ic=1,NC
    444           fluxtop(ic) = tablefluxtop(ic,j,isaisonflux)  *(1-factflux)
    445      .                + tablefluxtop(ic,j,isaisonflux+1)*  factflux
    446           fluxtop(ic) = fluxtop(ic)/((rinter(klev)-rinter(klev-1))*1e5)
    447         enddo
    448 
    449 c TEST: sortie de l'un des flux avec saveLs et factflux
    450 c       dans un fichier special pour tracer avec gnuplot
    451 c       if (j.eq.2) then
    452 c       open(unit=11,file="flux_80N.txt",status='old',position='append')
    453 c         write(11,*) ls_rad*180./RPI, factflux,
    454 c    .           ((rinter(llm)-rinter(llm-1))*1e5), fluxtop
    455 c         write(11,*) ls_rad*180./RPI, factflux,
    456 c    .     fluxtop(utilaer(3)+1)*((rinter(llm)-rinter(llm-1))*1e5),  !C2H2
    457 c    .     fluxtop(utilaer(6)+1)*((rinter(llm)-rinter(llm-1))*1e5)   !HCN
    458 c       close(11)
    459 c       endif
    460 c       if (j.eq.17) then
    461 c       open(unit=11,file="flux_eq.txt",status='old',position='append')
    462 c         write(11,*) ls_rad*180./RPI, factflux,
    463 c    .           ((rinter(llm)-rinter(llm-1))*1e5), fluxtop
    464 c         write(11,*) ls_rad*180./RPI, factflux,
    465 c    .     fluxtop(utilaer(3)+1)*((rinter(llm)-rinter(llm-1))*1e5),  !C2H2
    466 c    .     fluxtop(utilaer(6)+1)*((rinter(llm)-rinter(llm-1))*1e5)   !HCN
    467 c       close(11)
    468 c       endif
    469 c FIN TEST: sortie de l'un des flux avec ls et factflux
    470 
    471 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     382c extension pour klev+1 a NLEV avec tcp.ver
     383       ialt=1
     384       do l=klev+1,NLEV
     385           zalt=rmil(l)-RA/1000.
     386           do i=ialt,130
     387            if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
     388              ialt=i
     389            endif
     390           enddo
     391           factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
     392           press_c(l) = exp(  log(p1d(ialt))   * (1-factalt)
     393     &                      + log(p1d(ialt+1)) * factalt    )
     394           nb(l)      = exp(  log(ct1d(ialt))   * (1-factalt)
     395     &                      + log(ct1d(ialt+1)) * factalt    )
     396           temp_c(l)  = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
     397       enddo
     398               
     399c-----------------------------------------------------------------------
     400c
     401c   passage krpd => krate
     402c   ---------------------
     403c initialisation krate pour dissociations
     404
     405      if ((declin_c*10+267).lt.14.) then
     406          declinint = 0
     407          dec       = 0
     408      else
     409       if ((declin_c*10+267).gt.520.) then
     410          declinint = 14
     411          dec       = 534
     412       else
     413          declinint = 1
     414          dec       = 27
     415          do while( (declin_c*10+267).ge.real(dec+20) )
     416            dec       = dec+40
     417            declinint = declinint+1
     418          enddo
     419       endif
     420      endif
     421      if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
     422          factdec = ( declin_c - (dec-267)/10. ) / 4.
     423      else
     424          factdec = ( declin_c - (dec-267)/10. ) / 2.7
     425      endif
     426
     427      do l=1,NLEV
     428
     429c INTERPOL EN ALT POUR k (krpd tous les deux km)
     430        ialt    = int((rmil(l)-RA/1000.)/2.)+1
     431        factalt = (rmil(l)-RA/1000.)/2.-(ialt-1)
     432
     433        do i=1,ND+1
     434          krpddecm1 = krpd(declinint  ,i,ialt  ,klat(j)) * (1-factalt)
     435     &              + krpd(declinint  ,i,ialt+1,klat(j)) * factalt
     436          krpddec   = krpd(declinint+1,i,ialt  ,klat(j)) * (1-factalt)
     437     &              + krpd(declinint+1,i,ialt+1,klat(j)) * factalt
     438          krpddecp1 = krpd(declinint+2,i,ialt  ,klat(j)) * (1-factalt)
     439     &              + krpd(declinint+2,i,ialt+1,klat(j)) * factalt
     440
     441                    ! ND+1 correspond a la dissociation de N2 par les GCR
     442          if ( factdec.lt.0. ) then
     443        krate(l,i) = krpddecm1 * abs(factdec)
     444     &             + krpddec   * ( 1 + factdec)
     445          endif
     446          if ( factdec.gt.0. ) then
     447        krate(l,i) = krpddecp1 * factdec
     448     &             + krpddec   * ( 1 - factdec)
     449          endif
     450          if ( factdec.eq.0. ) krate(l,i) = krpddec
     451        enddo       
     452      enddo       
    472453
    473454c-----------------------------------------------------------------------
    … …  
    475456c   composition
    476457c   ------------
    477        
     458
    478459       do ic=1,NC
    479460        do l=1,klev
    480           cqy(l,ic) = qy_c(j,l,ic)
    481         enddo
    482        enddo
    483              
     461          cqy(l,ic)  = qy_c(j,l,ic)
     462          cqy0(l,ic) = cqy(l,ic)
     463        enddo
     464       enddo
     465
     466c lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a NLEV
     467
     468      WRITE(str2,'(i2.2)') klat(j)
     469      fich = "comp_"//str2//".dat"
     470      inquire (file=fich,exist=okfic)
     471      if (okfic) then
     472       open(11,file=fich,status='old')
     473c premiere ligne=declin
     474       read(11,'(A15)') ficdec
     475       write(curdec,'(E15.5)') declin_c
     476c si la declin est la meme: ce fichier a deja ete reecrit
     477c on lit la colonne t-1 au lieu de la colonne t
     478c (cas d une bande de latitude partagee par 2 procs)
     479       do ic=1,NC
     480        read(11,'(A10)') name
     481        if (name.ne.nomqy_c(ic)) then
     482          print*,"probleme lecture ",fich
     483          print*,name,nomqy_c(ic)
     484          stop
     485        endif
     486        if (ficdec.eq.curdec) then
     487          do l=klev+1,NLEV
     488            read(11,*) ficalt,cqy(l,ic),newq
     489          enddo
     490        else
     491          do l=klev+1,NLEV
     492            read(11,*) ficalt,oldq,cqy(l,ic)
     493          enddo
     494        endif
     495       enddo
     496       close(11)
     497      else       ! le fichier n'est pas la
     498       do ic=1,NC
     499        do l=klev+1,NLEV
     500          cqy(l,ic)=cqy(klev,ic)
     501        enddo
     502       enddo
     503      endif
     504      cqy0 = cqy
     505
    484506c-----------------------------------------------------------------------
    485507c
    … …  
    502524       
    503525       call gptitan(rinter,temp_c,nb,
    504      $              nomqy_c,cqy,fluxtop,
    505      $              declin_c,duree,(klat(j)-1),mass,
    506      $              botCH4,krpd,krate,reactif,
     526     $              nomqy_c,cqy,
     527     $              duree,(klat(j)-1),mass,md,
     528     $              kedd,botCH4,krate,reactif,
    507529     $              nom_prod,nom_perte,prod,perte,
    508530     $              aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,
    … …  
    514536       do ic=1,NC
    515537        do l=1,klev
    516           dqyc(j,l,ic) = (cqy(l,ic) - qy_c(j,l,ic))/dtchim  ! en /s
     538          dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim  ! en /s
    517539        enddo
    518540       enddo
    … …  
    535557
    536558       endif
     559
     560c-----------------------------------------------------------------------
     561c
     562c   sauvegarde compo sur NLEV
     563c   -------------------------
     564
     565c dans fichier compo_klat(j) (01 à 49)
     566       
     567      open(11,file=fich,status='replace')
     568c premiere ligne=declin
     569      write(11,'(E15.5)') declin_c
     570      do ic=1,NC
     571       write(11,'(A10)') nomqy_c(ic)
     572       do l=klev+1,NLEV
     573        write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-RA/1000.,
     574     .                                cqy0(l,ic),cqy(l,ic)
     575       enddo
     576      enddo
     577      close(11)
    537578
    538579c***********************************************************************

  • trunk/LMDZ.TITAN/libf/phytitan/common_mod.F90

    r1056 r1126  
    3434
    3535! ancien titan_for.h
    36       INTEGER, PARAMETER :: NLEV=llm,NC=44,ND=54,NR=377
     36      INTEGER, PARAMETER :: NLEV=llm+70,NC=44,ND=54,NR=377
    3737!$OMP THREADPRIVATE(NLEV,NC,ND,NR)
    3838!!!  doivent etre en accord avec titan.h
     39! pour l'UV (650 niveaux de 2 km)
     40      INTEGER, PARAMETER :: NLRT=650
    3941
    4042! ancien diagmuphy.h

  • trunk/LMDZ.TITAN/libf/phytitan/optci.F

    r1080 r1126  
    162162          XICLDI2(K)=TNI
    163163 420  CONTINUE
     164
     165c DEBUG
     166c       print*,"wnoi=",WNOI
    164167
    165168C
    … …  
    284287        TauGID(ig,:,:) = TAUGID_1pt(:,:)
    285288
     289c DEBUG
     290c     if(ig.eq.(ngrid/2+16)) then
     291c         print*,ig,'/',KLON,':'
     292c         print*,'TauHID_1',TAUHID(ig,1,:)
     293c         print*,'TauGID_1',TAUGID(ig,1,:)
     294c         print*,'TauHID_50',TAUHID(ig,50,:)
     295c         print*,'TauGID_50',TAUGID(ig,50,:)
     296c         print*,"DTAUI_1=",DTAUI(ig,1,:)
     297c         print*,"DTAUI_50=",DTAUI(ig,50,:)
     298c         print*,'cosBI_1',COSBI(ig,1,:)
     299c         print*,'cosBI_50',COSBI(ig,50,:)
     300c         print*,'WBARI_1',WBARI(ig,1,:)
     301c         print*,'WBARI_50',WBARI(ig,50,:)
     302c         stop
     303c     endif
     304
    286305c************************************************************************
    287306c************************************************************************

  • trunk/LMDZ.TITAN/libf/phytitan/optci_1pt_3.F

    r1058 r1126  
    282282c     CBAR=xv3(j,k)
    283283
     284c     print*, 'HERE, CIRS AEROSOLS'
     285c     call cirs_haze(PRESS(J),WNOI(K),TAEROS,TAEROSCAT,CBAR)
    284286
    285287c*********** EN TRAVAUX ***************************

  • trunk/LMDZ.TITAN/libf/phytitan/optcv.F

    r1081 r1126  
    148148c       close(1)
    149149
     150c DEBUG
     151c       print*,"wnov=",WNOV
     152
    150153      endif    ! fin initialisations premier appel
    151154
    … …  
    167170        else
    168171         if (ig.eq.1)  then
    169 c initialisation zqaer_1pt a partir d'une look-up table (uniforme en ig)
     172c initialisation zqaer_1pt a partir d une look-up table (uniforme en ig)
    170173c boucle sur nrad=10
    171174           open(10,file="qaer_eq_1d.dat")
    … …  
    220223        TauGVD(ig,:,:) = TAUGVD_1pt(:,:)
    221224
     225c DEBUG
     226c     if(ig.eq.(ngrid/2+16)) then
     227c         print*,ig,'/',KLON,':'
     228c         print*,'TauHVD_1',TAUHVD(ig,1,:)
     229c         print*,'TauGVD_1',TAUGVD(ig,1,:)
     230c         print*,'TauHVD_50',TAUHVD(ig,50,:)
     231c         print*,'TauGVD_50',TAUGVD(ig,50,:)
     232c     stop
     233c     endif
     234
    222235 101  CONTINUE
    223236

  • trunk/LMDZ.TITAN/libf/phytitan/optcv_1pt_3.F

    r1058 r1126  
    231231
    232232       IF (TAEROSCAT.le.0.) CBAR=0.
     233
     234c     print*, 'HERE, CIRS AEROSOLS'
     235c     call cirs_haze(PRESS(J),WNOV(K),TAEROS,TAEROSCAT,CBAR)
    233236
    2342371699  FORMAT(a3,2I3,3(ES15.7,1X))

  • trunk/LMDZ.TITAN/libf/phytitan/physiq.F

    r1120 r1126  
    419419         itap    = 0
    420420         itaprad = 0
    421          itapchim    = 0
     421         itapchim    = 1
    422422
    423423c init rnuabar
    … …  
    11261126         tr_seri(:,:,1:nmicro) = tr_seri(:,:,1:nmicro)
    11271127     .                        + d_tr_mph(:,:,1:nmicro)*dtime
    1128 
     1128c        call WriteField_phy('physiq_d_tr_mph01',
     1129c    .                          d_tr_mph(:,:,1),klev)
     1130c        call WriteField_phy('physiq_d_tr_mph10',
     1131c    .                          d_tr_mph(:,:,10),klev)
     1132        endif
    11291133c       PAS ELEGANT mais je n'ai pas trouve d'autres solutions :
    11301134c       Il semblerait qu'il y ait un probleme lorsque les tendances de traceurs
    … …  
    11411145        ENDDO
    11421146
    1143         endif
    1144 
    11451147c condensation:
    11461148c       NE PAS OUBLIER LA CONDENSATION DES NUAGES !!!!
    … …  
    11491151     .                         + d_tr_mph(:,:,nmicro+1:nqmax)*dtime
    11501152        endif
     1153
     1154c chimie:
    11511155        if ((chimi).and.(nqmax.gt.nmicro)) then
    1152 c chimie:
    11531156         tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_kim(:,:,:)*dtime
    11541157        endif
    11551158
    1156       endif
     1159      endif  !iflag_trac=1
    11571160
    11581161c       ch4=0.

  • trunk/LMDZ.TITAN/libf/phytitan/phytrac.F

    r1072 r1126  
    459459c ------------
    460460       CALL calchim(klon,nqmax-nmicro,ychim,nomqy,pdecli,lonsol,dtkim,
    461      .              ztemp,zplay,zplev,
     461     .              ztemp,zplay,zplev,zzlay,zzlev,
    462462     .              pdyfi)   
    463463c ychim ne doit pas etre modifie, pdyfi en /s

  • trunk/LMDZ.TITAN/libf/phytitan/profile.F

    r904 r1126  
    9494       b2 =       3183.
    9595       c2 =       4737.
     96c pres est en Pa => conversion car expression veut p en hPa
    9697       DO il=1,nlev
    97          temp(il)=a1*exp(-((pres(il)-b1)/c1)**2.)
    98      .          + a2*exp(-((pres(il)-b2)/c2)**2.)
     98         temp(il)=a1*exp(-((pres(il)/100.-b1)/c1)**2.)
     99     .          + a2*exp(-((pres(il)/100.-b2)/c2)**2.)
    99100       ENDDO
    100101       zkm(1)  = 0.0
Note: See TracChangeset for help on using the changeset viewer.