Changeset 1908


Ignore:
Timestamp:
Mar 8, 2018, 5:16:35 PM (7 years ago)
Author:
jvatant
Message:

Making chemistry handling more flexible - Major and Final Step (hopefully) !
After preliminary commits r1871-86-87-91-94-98, here comes the major update of the interface
with photochemical module + fix how tendencies for chem and mufi tracers are managed in physiq_mod !
+ Major modifs in :
++ calchim.F90 to comply with flexible resolution, parallelism, upper pressure grid ...
++ physiq_mod.F90 where there was a bug on the update of the tracers and their tendencies for calchim

and calmufi ( we actually were sending non-updated fields to these processes )
We also now put the same tendency on all longitudes within a lat band and not
relative tendency if 2D chemistry ( and we set to zero if ever negs are created )

+ Also modifs to have chemistry in 1D in rcm1d ( and moved gr_kim_vervack in phytitan to be accessible for 1d )
+ In chemistry added a check.c to verify coherence of sizes between C and Fortran
--JVO

Location:
trunk/LMDZ.TITAN
Files:
2 added
1 deleted
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/README

    r1904 r1908  
    13991399+ Upper chemistry fields ( and chemistry tracers if needed ) in newstart are
    14001400 regrided in vert_regrid_kim.F90
     1401+ Interface with chemistry routines itself will be done in a later commit (r1908)
    14011402
    14021403== 13/01/18 == JVO
    … …  
    14231424== 25/02/18 == JVO
    14241425Fix two bugs, one related to zonal means and one to rat_mmol
    1425 + Anyway zonal means apparently lead to crashes in microphysics and chemistry.
    1426 + Chemistry is still on the bench, but for microphysics, please use full 3D !!
     1426+ Zonal means may lead to crashes in microphysics, please use full 3D, which is very quick !!
     1427
     1428== 08/03/18 == JVO
     1429Making chemistry handling more flexible - Major and Final Step (hopefully) !
     1430After preliminary commits r1871-86-87-91-94-98, here comes the major update of the interface
     1431with photochemical module + fix how tendencies for chem and mufi tracers are managed in physiq_mod !
     1432+ Major modifs in :
     1433++ calchim.F90 to comply with flexible resolution, parallelism, upper pressure grid ...
     1434++ physiq_mod.F90 where there was a bug on the update of the tracers and their tendencies for calchim
     1435   and calmufi ( we actually were sending non-updated fields to these processes )
     1436   We also now put the same tendency on all longitudes within a lat band and not
     1437   relative tendency if 2D chemistry ( and we set to zero if ever negs are created )
     1438+ Also modifs to have chemistry in 1D in rcm1d ( and moved gr_kim_vervack in phytitan to be accessible for 1d )
     1439+ In chemistry added a check.c to verify coherence of sizes between C and Fortran

  • trunk/LMDZ.TITAN/libf/chimtitan/gptitan.c

    r1126 r1908  
    11/* gptitan: photochimie */
    22/* GCCM */
    3 
    4 /* tout est passe en simple precision */
    5 /* sauf pour l'inversion de la matrice */
    63
    74/* nitriles et hydrocarbures separes pour l'inversion */
    … …  
    1411     double *RA, double *TEMP, double *NB,
    1512     char CORPS[][10], double Y[][NLEV],
    16      double *FIN, int *LAT, double *MASS, double MD[][NLEV],
     13     double *FIN, double *LAT, double *MASS, double MD[][NLEV],
    1714     double *KEDD, double *botCH4, double KRATE[][NLEV],
    1815     int reactif[][5], int *nom_prod, int *nom_perte,
    … …  
    5350
    5451/* DEBUG */
    55       printf("CHIMIE: lat=%d\n",(*LAT)+1);
     52      printf("CHIMIE: lat=%g\n",(*LAT));
    5653/**/
    5754
    … …  
    6562   strcat( outlog, ".log" );
    6663   out = fopen( outlog, "w" );
    67    fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
     64   fprintf(out,"CHIMIE: lat=%g\n",(*LAT));
    6865   fclose( out );
    6966
    … …  
    7976/*  DEBUG
    8077            out = fopen( outlog, "a" );
    81    fprintf(out,"CHIMIE: lat=%d\n",(*LAT)+1);
     78   fprintf(out,"CHIMIE: lat=%g\n",(*LAT));
    8279            fclose( out );
    8380*/
    … …  
    236233             }
    237234/* DEBUG
    238 printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
     235printf("AERPROD : LAT = %g - J = %d\n",(*LAT),j);
    239236if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
    240237      printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
    … …  
    286283
    287284/* DEBUG
    288 printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
     285printf("HTOH2 : LAT = %g - J = %d\n",(*LAT),j);
    289286if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
    290287printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
    … …  
    511508                  {
    512509                     out = fopen( outlog, "a" );
    513                      fprintf( out, "Lat no %d;", (*LAT)+1);
     510                     fprintf( out, "Latitude %g;", (*LAT));
    514511                     fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
    515512                     fclose( out );
    … …  
    666663             }
    667664/* DEBUG
    668 printf("AERPROD : LAT = %d - J = %d\n",(*LAT),j);
     665printf("AERPROD : LAT = %g - J = %d\n",(*LAT),j);
    669666if(fabs(dyc2h2*NB[j])>fabs(fp[utilaer[2]]/10.))
    670667      printf("fp(%s) =%e; dyc2h2 =%e\n",corps[utilaer[2]],
    … …  
    716713
    717714/* DEBUG
    718 printf("HTOH2 : LAT = %d - J = %d\n",(*LAT),j);
     715printf("HTOH2 : LAT = %g - J = %d\n",(*LAT),j);
    719716if(fabs(dyh*NB[j])>fabs(fp[utilaer[0]]/10.))
    720717printf("fp(%s) = %e; dyh  = %e\n",corps[utilaer[0]],fp[utilaer[0]],dyh*NB[j]);
    … …  
    810807                  {
    811808                     out = fopen( outlog, "a" );
    812                      fprintf( out, "Lat no %d; declin:%e;", (*LAT)+1, (*DECLIN) );
     809                     fprintf( out, "Latitude: %g; declin:%e;", (*LAT), (*DECLIN) );
    813810                     fprintf(out, " alt:%e; %s %e %e ; %e %e\n",(RA[j]-R0),corps[i],ym1[i],Y[i][j],time,delta);
    814811                     fclose( out );

  • trunk/LMDZ.TITAN/libf/chimtitan/titan.h

    r1126 r1908  
    77
    88#define R0    (double)(2575.0) /* Titan's radius */
    9 #define NLEV  (int)(125)  /* Nbre de niv verticaux - =llm+70 dans common_mod */
    10 #define NLD   (int)(40)   /* Nbre de niv verticaux faits sans diff */
     9#define NLEV  (int)(133)  /* Nbre de niv verticaux - =llm+70 dans common_mod -> Need to be coherent with the vertical grid used !! */
     10#define NLD   (int)(40)   /* Nbre de niv verticaux faits sans diff -> Need to be coherent with the vertical grid used !! */
    1111#define NLRT  (int)(650)  /* Nbre de niv verticaux dans table fmoy - aussi dans common_mod */
    1212

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

    r1903 r1908  
    1       SUBROUTINE calchim(ngrid,qy_c,declin_rad,ls_rad, &
    2                          dtchim,ctemp,cplay,cplev,czlay,czlev,dqyc)
    3      
    4 !--------------------------------------------------------------------
    5 ! 
    6 !     Introduction d une routine chimique
    7 !
    8 !     Auteurs: S. Lebonnois,  01/2000 | 09/2003
    9 !            adaptation pour Titan 3D: 02/2009
    10 !            adaptation pour // : 04/2013
    11 !            extension chimie jusqu a 1300km : 10/2013
    12 !     
    13 !              J. Vatant d'Ollone, 02/2017
    14 !               + adaptation pour le nouveau titan issu du generic
    15 !
    16 !---------------------------------------------------------------------
    17 !
    18      
    19       use comchem_h
    20       use dimphy
    21       use datafile_mod, only: datadir
    22       use comcstfi_mod, only: rad, pi, kbol
    23       use moyzon_mod, only: tmoy,playmoy,zlaymoy,zlevmoy,klat
    24       use mod_grid_phy_lmdz, only: nbp_lat
    25      
    26       implicit none
     1SUBROUTINE calchim(ngrid,qy_c,declin,dtchim,            &
     2     ctemp,cpphi,cplay,cplev,czlay,czlev,dqyc)
     3
     4  !---------------------------------------------------------------------------------------------------------
     5  ! 
     6  ! Purpose : Interface subroutine to photochemical C model for Titan GCM.
     7  ! -------
     8  !           The subroutine computes the chemical processes for a single vertical column.
     9  !
     10  ! - Only tendencies are returned.
     11  ! - With moyzon_ch=.true. and input vectors zonally averaged
     12  !   the calculation is done only once per lat. band
     13  !
     14  ! Authors: + S. Lebonnois : 01/2000 | 09/2003
     15  ! -------                  adaptation for Titan 3D : 02/2009
     16  !                          adaptation for // : 04/2013
     17  !                          extension chemistry up to 1300km : 10/2013
     18  !     
     19  !          + J. Vatant d'Ollone
     20  !               + 02/17 - adaptation for the issu du generic
     21  !               + 01/18 - 03/18 - Major transformations :
     22  !                   - Upper chemistry fields are now stored in startfi
     23  !                     and defined on a pressure grid from Vervack profile
     24  !                   - Altitudes above GCM top are computed using correct g(z)
     25  !                     to be coherent with obs and with dz=10km for UV processes
     26  !                     but this should be even better if all physiq "under" was doing this !
     27  !                   - These modifs enables to run chemistry with others resolution than 32x48x55 !
     28  !                   - Only the actinic fluxes are still read in a 49-lat input but interp. on lat grid
     29  !                   - Chemistry can still be done in 2D
     30  !                      -> Calcul. once per band lat and put same tendency in all longi.
     31  !                         Check for negs in physiq_mod.
     32  !                      -> If procs sharing a lat band, no problem, the calcul will just be done twice.
     33  !                      -> Will not work with Dynamico, where the chemistry will have to be done in 3D.
     34  !                          ( and there'll be work to do to get rid of averaged fields )
     35  !
     36  ! + STILL TO DO : Replug the interaction with haze (cf titan.old) -> to see with JB.
     37  !---------------------------------------------------------------------------------------------------------
     38
     39  ! --------------------------------------------------------------------
     40  ! Structure :
     41  ! -----------
     42  !   0.  Declarations
     43  !   I.  Init and firstcall
     44  !         1. Read and store Vervack profile
     45  !         2. Compute planetar averaged atm. properties
     46  !         3. Init compound caracteristics
     47  !         4. Init photodissociations rates from actinic fluxes
     48  !         5. Init chemical reactions
     49  !         6. Init eddy diffusion coeff
     50  !   II. Loop on latitudes/grid-points
     51  !         0. Check on 2D chemistry
     52  !           1. Compute atm. properties at grid points
     53  !           2. Interpolate photodissociation rates at lat,alt,dec
     54  !           3. Read composition
     55  !           4. Call main solver gptitan C routine
     56  !           5. Calculate output tendencies on advected tracers
     57  !           6. Update upper chemistry fields
     58  !         0bis. If 2D chemsitry, don't recalculate if needed
     59  ! -----------------------------------------------------------------
     60
     61  USE comchem_h
     62  USE dimphy
     63  USE datafile_mod, ONLY: datadir
     64  USE comcstfi_mod, ONLY: g, rad, pi, r, kbol
     65  USE geometry_mod, ONLY: latitude
     66  USE logic_mod, ONLY: moyzon_ch
     67  USE moyzon_mod, ONLY: tmoy, playmoy, phimoy, zlaymoy, zlevmoy
     68
     69  IMPLICIT NONE
    2770
    2871! ------------------------------------------
    2972! *********** 0. Declarations *************
    3073! ------------------------------------------
    31      
    32 !    Arguments
    33 !    ---------
    34 
    35       INTEGER      ngrid                   ! nb of horiz points
    36       REAL*8         qy_c(ngrid,klev,nkim)     ! Especes chimiques apres adv.+diss.
    37       REAL*8         declin_rad,ls_rad      ! declinaison et long solaire en radians
    38       REAL*8         dtchim                 ! pas de temps chimie
    39       REAL*8         ctemp(ngrid,klev)       ! Temperature
    40       REAL*8         cplay(ngrid,klev)       ! pression (Pa)
    41       REAL*8         cplev(ngrid,klev+1)     ! pression intercouches (Pa)
    42       REAL*8         czlay(ngrid,klev)       ! altitude (m)
    43       REAL*8         czlev(ngrid,klev+1)     ! altitude intercouches (m)
    44      
    45       REAL*8         dqyc(ngrid,klev,nkim)     ! Tendances especes chimiques
    46      
    47      
    48 !    Local variables :
    49 !    -----------------
    50 
    51       integer i,j,l,ic,jm1
    52 
    53 ! variables envoyees dans la chimie: double precision
    54 
    55       REAL*8  temp_c(nlaykim_tot)
    56       REAL*8  press_c(nlaykim_tot)   ! T,p(mbar) a 1 lat donnee
    57       REAL*8  cqy(nlaykim_tot,nkim)    ! frac mol qui seront modifiees
    58       REAL*8  cqy0(nlaykim_tot,nkim)    ! frac mol avant modif
    59       REAL*8  surfhaze(nlaykim_tot)
    60       REAL*8  cprodaer(nlaykim_tot,4),cmaer(nlaykim_tot,4)
    61       REAL*8  ccsn(nlaykim_tot,4),ccsh(nlaykim_tot,4)
    62 ! rmil: milieu de couche, grille pour K,D,p,ct,T,y
    63 ! rinter: intercouche (grille RA dans la chimie)
    64       REAL*8  rmil(nlaykim_tot),rinter(nlaykim_tot),nb(nlaykim_tot)
    65       REAL*8,save,allocatable :: kedd(:)
    66 
    67       character str1*1,str2*2
    68      
    69 !     variables locales initialisees au premier appel
    70 
    71       LOGICAL firstcall
    72       DATA    firstcall/.true./
    73       SAVE    firstcall
    74      
    75       integer dec,declinint,ialt
    76       REAL*8  declin_c                       ! declinaison en degres
    77       real*8  factalt,factdec,krpddec,krpddecp1,krpddecm1
    78       real*8  duree
    79       REAL*8,save :: mass(nkim)
    80       REAL*8,save,allocatable :: md(:,:)
    81       REAL*8,save :: botCH4
    82       DATA  botCH4/0.05/
    83       real*8,save ::  r1d(131),ct1d(131),p1d(131),t1d(131) ! lecture tcp.ver
    84       REAL*8,save,allocatable :: krpd(:,:,:,:),krate(:,:)
    85       integer,save :: reactif(5,nr_kim),nom_prod(nkim),nom_perte(nkim)
    86       integer,save :: prod(200,nkim),perte(2,200,nkim)
    87       character  fich*7,ficdec*15,curdec*15,name*10
    88       real*8  ficalt,oldq,newq,zalt
    89       logical okfic
    90 
    91 !     Dummy parameters without microphysics
    92 !      + Just here to keep the whole stuff running without modif C sources
    93 !     ---------------------------------------------------------------------
    94 
    95       INTEGER  :: utilaer(16)
    96       INTEGER  :: aerprod = 0
    97       INTEGER  :: htoh2   = 0
    98 
    99 !-----------------------------------------------------------------------
    100 !***********************************************************************
    101 !
    102 !    Initialisations :
    103 !    ----------------
    104 
    105       duree = dtchim  ! passage en real*4 pour appel a routines C
    106 
    107       IF (firstcall) THEN
    108            
    109         print*,'CHIMIE, premier appel'
    110        
    111 ! ************************************
    112 ! Au premier appel, initialisation de toutes les variables
    113 ! pour les routines de la chimie.
    114 ! ************************************
    115 
    116         allocate(kedd(nlaykim_tot))
    117 
    118         allocate(krpd(15,nd_kim+1,nlrt_kim,nbp_lat),krate(nlaykim_tot,nr_kim),md(nlaykim_tot,nkim))
    119 
    120 
    121 ! calcul de temp_c, densites et press_c en moyenne planetaire :
    122 ! --------------------------------------------------------------
    123 
    124         print*,'pression, densites et temp (init chimie):'
    125         print*,'level, press_c, nb, temp_c'
     74
     75  ! Arguments
     76  ! ---------
     77
     78  INTEGER, INTENT(IN)                              :: ngrid       ! Number of atmospheric columns.
     79  REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(IN)   :: qy_c        ! Chemical species on GCM layers after adv.+diss. (mol/mol).
     80  REAL*8, INTENT(IN)                               :: declin      ! Solar declination (rad).
     81  REAL*8, INTENT(IN)                               :: dtchim      ! Chemistry timsetep (s).
     82  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: ctemp       ! Mid-layer temperature (K).
     83  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: cpphi       ! Mid-layer geopotential (m2.s-2).
     84  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: cplay       ! Mid-layer pressure (Pa).
     85  REAL*8, DIMENSION(ngrid,klev+1),    INTENT(IN)   :: cplev       ! Inter-layer pressure (Pa).
     86  REAL*8, DIMENSION(ngrid,klev),      INTENT(IN)   :: czlay       ! Mid-layer altitude (m).
     87  REAL*8, DIMENSION(ngrid,klev+1),    INTENT(IN)   :: czlev       ! Inter-layer altitude (m).
     88
     89  REAL*8, DIMENSION(ngrid,klev,nkim), INTENT(OUT)  :: dqyc        ! Chemical species tendencies on GCM layers (mol/mol/s).
     90
     91  ! Local variables :
     92  ! -----------------
     93
     94  INTEGER :: i , l, ic, ig, igm1
     95
     96  INTEGER :: dec, idec, ipres, ialt, klat
     97
     98  REAL*8  :: declin_c  ! Declination (deg).
     99  REAL*8  :: factp, factalt, factdec, factlat, krpddec, krpddecp1, krpddecm1
     100  REAL*8  :: temp1, logp
     101
     102  ! Variables sent into chemistry module (must be in double precision)
     103  ! ------------------------------------------------------------------
     104
     105  REAL*8, DIMENSION(nlaykim_tot) :: temp_c  ! Temperature (K).
     106  REAL*8, DIMENSION(nlaykim_tot) :: press_c ! Pressure (Pa).
     107  REAL*8, DIMENSION(nlaykim_tot) :: phi_c   ! Geopotential (m2.s-2) - actually not sent in chem. module but used to compute alts.
     108  REAL*8, DIMENSION(nlaykim_tot) :: nb      ! Density (cm-3).
     109
     110  REAL*8, DIMENSION(nlaykim_tot,nkim) :: cqy   ! Chemical species in whole column (mol/mol) sent to chem. module.
     111  REAL*8, DIMENSION(nlaykim_tot,nkim) :: cqy0  !     "      "     "   "      "        "     before modifs.
     112
     113  REAL*8  :: surfhaze(nlaykim_tot)
     114  REAL*8  :: cprodaer(nlaykim_tot,4), cmaer(nlaykim_tot,4)
     115  REAL*8  :: ccsn(nlaykim_tot,4), ccsh(nlaykim_tot,4)
     116
     117  REAL*8, DIMENSION(nlaykim_tot) :: rmil   ! Mid-layer distance (km) to planetographic center.
     118  REAL*8, DIMENSION(nlaykim_tot) :: rinter ! Inter-layer distance (km) to planetographic center (RA grid in chem. module).
     119  ! NB : rinter is on nlaykim_tot too, we don't care of the uppermost layer upper boundary altitude.
     120
     121  ! Saved variables initialized at firstcall
     122  ! ----------------------------------------
     123
     124  LOGICAL, SAVE :: firstcall = .TRUE.
     125!$OMP THREADPRIVATE(firstcall)
     126
     127  REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: kedd ! Eddy mixing coefficient for upper chemistry (cm^2.s-1)
     128!$OMP THREADPRIVATE(kedd)
     129
     130  REAL*8, DIMENSION(:),   ALLOCATABLE, SAVE :: mass ! Molar mass of the compunds (g.mol-1)
     131  REAL*8, DIMENSION(:,:), ALLOCATABLE, SAVE :: md
     132!$OMP THREADPRIVATE(mass,md)
     133
     134  REAL*8, DIMENSION(:), ALLOCATABLE, SAVE :: r1d, ct1d, p1d, t1d ! Vervack profile
     135  ! JVO 18 : No threadprivate for those as they'll be read in tcp.ver by master
     136
     137  REAL*8, DIMENSION(:,:,:,:), ALLOCATABLE, SAVE :: krpd  ! Photodissociations rate table
     138  REAL*8, DIMENSION(:,:)    , ALLOCATABLE, SAVE :: krate ! Reactions rate ( photo + chem )
     139!$OMP THREADPRIVATE(krpd,krate)
     140
     141  INTEGER, DIMENSION(:),     ALLOCATABLE, SAVE :: nom_prod, nom_perte
     142  INTEGER, DIMENSION(:,:),   ALLOCATABLE, SAVE :: reactif
     143  INTEGER, DIMENSION(:,:),   ALLOCATABLE, SAVE :: prod
     144  INTEGER, DIMENSION(:,:,:), ALLOCATABLE, SAVE :: perte
     145!$OMP THREADPRIVATE(nom_prod,nom_perte,reactif,prod,perte)
     146
     147  ! TEMPORARY : Dummy parameters without microphysics
     148  ! Here to keep the whole stuff running without modif chem. module
     149  ! ---------------------------------------------------------------
     150
     151  INTEGER  :: utilaer(16)
     152  INTEGER  :: aerprod = 0
     153  INTEGER  :: htoh2   = 0
     154
     155  ! -----------------------------------------------------------------------
     156  ! ***************** I. Initialisations and Firstcall ********************
     157  ! -----------------------------------------------------------------------
     158
     159  IF (firstcall) THEN
     160
     161     PRINT*, 'CHIMIE, premier appel'
     162
     163     call check(nlaykim_tot,klev-15,nlrt_kim,nkim)
     164
     165     ALLOCATE(r1d(131))
     166     ALLOCATE(ct1d(131))
     167     ALLOCATE(p1d(131))
     168     ALLOCATE(t1d(131))
     169
     170     ALLOCATE(md(nlaykim_tot,nkim))
     171     ALLOCATE(mass(nkim))
     172     
     173     ALLOCATE(kedd(nlaykim_tot))
     174     
     175     ALLOCATE(krate(nlaykim_tot,nr_kim))
     176     ALLOCATE(krpd(15,nd_kim+1,nlrt_kim,nlat_actfluxes))
     177
     178     ALLOCATE(nom_prod(nkim))
     179     ALLOCATE(nom_perte(nkim))
     180     ALLOCATE(reactif(5,nr_kim))
     181     ALLOCATE(prod(200,nkim))
     182     ALLOCATE(perte(2,200,nkim))
     183     
     184     ! 1. Read Vervack profile "tcp.ver", once for all
     185     ! -----------------------------------------------
     186
     187!$OMP MASTER
     188     OPEN(11,FILE=TRIM(datadir)//'/tcp.ver',STATUS='old')
     189     READ(11,*)
     190     DO i=1,131
     191        READ(11,*) r1d(i), t1d(i), ct1d(i), p1d(i)
     192        ! For debug :
     193        ! -----------
     194        ! PRINT*, "tcp.ver", r1d(i), t1d(i), ct1d(i), p1d(i)
     195     ENDDO
     196     CLOSE(11)
     197!$OMP END MASTER
     198!$OMP BARRIER
     199
     200     ! 2. Calculation of temp_c, densities and altitudes in planetary average
     201     ! ----------------------------------------------------------------------
     202     
     203     ! JVO18 : altitudes calculated in firstcall are just for diag, as I set kedd in pressure grid
     204
     205     ! a. For GCM layers we just copy-paste
     206
     207!$OMP MASTER
     208     PRINT*,'Init chemistry : pressure, density, temperature ... :'
     209     PRINT*,'level, rmil (km), press_c (mbar), nb (cm-3), temp_c(K)'
     210!$OMP END MASTER
     211!$OMP BARRIER
     212     
     213     IF (ngrid.NE.1) THEN
     214       DO l=1,klev
     215          rinter(l)  = (zlevmoy(l)+rad)/1000.0              ! km
     216          rmil(l)    = (zlaymoy(l)+rad)/1000.0              ! km
     217          temp_c(l)  = tmoy(l)                              ! K
     218          phi_c(l)   = phimoy(l)                            ! m2.s-2
     219          press_c(l) = playmoy(l)/100.                      ! mbar
     220          nb(l)      = 1.e-4*press_c(l) / (kbol*temp_c(l))  ! cm-3
     221!$OMP MASTER
     222          PRINT*, l, rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
     223!$OMP END MASTER
     224!$OMP BARRIER
     225       ENDDO
     226       rinter(klev+1)=(zlevmoy(klev+1)+rad)/1000.
     227     ELSE
     228       DO l=1,klev
     229          rinter(l)  = (czlev(1,l)+rad)/1000.0                ! km
     230          rmil(l)    = (czlay(1,l)+rad)/1000.0                ! km
     231          temp_c(l)  = ctemp(1,l)                             ! K
     232          phi_c(l)   = cpphi(1,l)                             ! m2.s-2
     233          press_c(l) = cplay(1,l)/100.                        ! mbar
     234          nb(l)      = 1.e-4*press_c(l) / (kbol*temp_c(l))  ! cm-3
     235          PRINT*, l, rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
     236       ENDDO
     237       rinter(klev+1)=(czlev(1,klev+1)+rad)/1000.
     238     ENDIF ! if ngrid.ne.1
     239
     240     ! b. Extension in upper atmosphere with Vervack profile
     241     !    NB : Maybe the transition klev/klev+1 is harsh ...     
     242
     243     ipres=1
     244     DO l=klev+1,nlaykim_tot
     245        press_c(l) = preskim(l-klev) / 100.0
     246        DO i=ipres,130
     247           IF ( (press_c(l).LE.p1d(i)) .AND. (press_c(l).GT.p1d(i+1)) ) THEN
     248              ipres=i
     249           ENDIF
     250        ENDDO
     251        factp = (press_c(l)-p1d(ipres)) / (p1d(ipres+1)-p1d(ipres))
     252
     253        nb(l)      = exp( log(ct1d(ipres))*(1-factp) + log(ct1d(ipres+1))* factp )
     254        temp_c(l)  = t1d(ipres)*(1-factp) + t1d(ipres+1)*factp
     255     ENDDO
     256
     257     ! We build altitude with hydrostatic equilibrium on preskim grid with Vervack profile
     258     ! ( keeping in mind that preskim is built based on Vervack profile with dz=10km )
     259
     260     DO l=klev+1,nlaykim_tot
     261
     262        ! Compute geopotential on the upper grid, to have correct altitude
     263        ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     264        ! NB : Upper chemistry model is based on observational profiles and needs correct altitudes !!
     265        ! Altitudes can (and must) be calculated with correct g even if it's not the case in the GCM under
     266
     267        temp1 = 0.5*(temp_c(l-1)+temp_c(l)) ! interlayer temp
     268        phi_c(l) = phi_c(l-1) + r*temp1*log(press_c(l-1)/press_c(l)) ! Geopotential assuming hydrostatic equilibrium
     269     
     270        rmil(l)  =  ( g*rad*rad / (g*rad - phi_c(l)) ) / 1000.0 ! z(phi) with g varying with altitude
     271
     272!$OMP MASTER
     273        PRINT*, l , rmil(l)-rad/1000.0, press_c(l), nb(l), temp_c(l), phi_c(l)
     274!$OMP END MASTER
     275!$OMP BARRIER
     276     ENDDO
     277
     278     DO l=klev+2,nlaykim_tot
     279        rinter(l) = 0.5*(rmil(l-1) + rmil(l))
     280     ENDDO
     281
     282     ! 3. Compounds caracteristics
     283     ! ---------------------------
     284     mass(:) = 0.0
     285     call comp(nomqy_c,nb,temp_c,mass,md)
     286     PRINT*,'           Mass'
     287     DO ic=1,nkim
     288        PRINT*, nomqy_c(ic), mass(ic)
     289     ENDDO
     290
     291     ! 4. Photodissociation rates
     292     ! --------------------------
     293     call disso(krpd,nlat_actfluxes)
     294
     295     ! 5. Init. chemical reactions with planetary average T prof.
     296     ! ----------------------------------------------------------
     297
     298     !  NB : Chemical reactions rate are assumed to be constant within the T range of Titan's atm
     299     !  so we fill their krate once for all but krate for photodiss will be filled at each timestep
     300     
     301     call chimie(nomqy_c,nb,temp_c,krate,reactif, &
     302          nom_perte,nom_prod,perte,prod)
     303
     304     ! 6. Eddy mixing coefficients (constant with time and space)
     305     ! ----------------------------------------------------------
     306
     307     kedd(:) = 1.e3 ! Default value =/= zero
     308
     309     ! NB : Eddy coeffs from 1D models (e.g. Lavvas et al 08) in altitude but they're rather linked to pressure
     310     !      Below GCM top we have dynamic mixing and for levs < nld=klev-15 the chem. solver ignores diffusion
     311
     312     ! First calculate kedd for upper chemistry layers
     313     DO l=klev+1,nlaykim_tot
     314        logp=-log10(press_c(l))
     315     ! 1E6 at 300 km ~ 10-1 mbar
     316        IF     ( logp.ge.1.0 .and. logp.le.4.0 ) THEN
     317              kedd(l) = 10.**(6.0+1.1*(logp-1.0))
     318     ! 2E7 at 600 km ~ 10-4 mbar
     319        ELSEIF ( logp.gt.4.0 .and. logp.le.6.0 ) THEN
     320              kedd(l) = 10.**(7.3+0.35*(logp-4.0))
     321     ! 1E8 above 900 km ~ 10-6 mbar
     322        ELSEIF ( logp.gt.6.0                   ) THEN
     323             kedd(l) = 1.e8
     324        ENDIF
     325     ENDDO
     326     
     327     ! Then adjust 15 last layers profile fading to default value depending on kedd(ptop)
     328     DO l=klev-15,klev
     329        temp1   = ( log10(press_c(l)/press_c(klev-15)) ) / ( log10(press_c(klev+1)/press_c(klev-15)) )
     330        kedd(l) = 10.**( 3.0 + log10(kedd(klev+1)/1.e3) * temp1 )
     331     ENDDO
     332     
     333     firstcall = .FALSE.
     334  ENDIF  ! firstcall
     335 
     336  declin_c = declin*180./pi
     337
     338  ! -----------------------------------------------------------------------
     339  ! *********************** II. Loop on latitudes *************************
     340  ! -----------------------------------------------------------------------
     341 
     342  DO ig=1,ngrid
     343
     344    IF (ig.eq.1) THEN
     345        igm1=1
     346     ELSE
     347        igm1=ig-1
     348     ENDIF
     349
     350     ! If 2D chemistry, trick to do the calculation only once per latitude band within the chunk
     351     ! NB1 : Will be obsolete with DYNAMICO, the chemistry will necessarly be 3D
     352     ! NB2 : Test of same latitude with dlat=0.5 : I think that if you run sims better than 1/10th degree then
     353     ! either it's with Dynamico and doesn't apply OR it is more than enough in terms of "preco / calc time" !
     354     ! -------------------------------------------------------------------------------------------------------
     355
     356     IF ( ( moyzon_ch .AND. ( ig.EQ.1 .OR. (ABS(latitude(ig)-latitude(igm1)).GT.0.1*pi/180.)) ) .OR. (.NOT. moyzon_ch) ) THEN
     357
     358        ! 1. Compute altitude for the grid point with hydrostat. equilib.
     359        ! ---------------------------------------------------------------
     360
     361        ! a. For GCM layers we just copy-paste
     362
    126363        DO l=1,klev
    127          rinter(l)  = (zlevmoy(l)+rad)/1000.
    128          rmil(l)    = (zlaymoy(l)+rad)/1000.
    129 !     temp_c (K):
    130          temp_c(l)  = tmoy(l)
    131 !     press_c (mbar):
    132          press_c(l) = playmoy(l)/100.
    133 !     nb (cm-3):
    134          nb(l) = 1.e-4*press_c(l) / (kbol*temp_c(l))
    135          print*,l,rmil(l)-rad/1000.,press_c(l),nb(l),temp_c(l)
    136         ENDDO
    137         rinter(klev+1)=(zlevmoy(klev+1)+rad)/1000.
    138 
    139 ! au-dessus du GCM, dr regulier et rinter(nlaykim_tot)=1290+2575 km.
    140        do l=klev+2,nlaykim_tot
    141          rinter(l) = rinter(klev+1) &
    142                + (l-klev-1)*(3865.-rinter(klev+1))/(nlaykim_tot-klev-1)
    143          rmil(l-1) = (rinter(l-1)+rinter(l))/2.
    144        enddo
    145        rmil(nlaykim_tot) = rinter(nlaykim_tot)+(rinter(nlaykim_tot)-rinter(nlaykim_tot-1))/2.
    146 ! rmil et rinter ne servent que pour la diffusion (> plafond-100km) donc ok en moyenne planetaire
    147 
    148 ! lecture de tcp.ver, une seule fois
    149 ! remplissage r1d,t1d,ct1d,p1d
    150         open(11,file=TRIM(datadir)//'/tcp.ver',status='old')
    151         read(11,*)
    152         do i=1,131
    153           read(11,*) r1d(i),t1d(i),ct1d(i),p1d(i)
    154           !print*,"TCP.VER ",r1d(i),t1d(i),ct1d(i),p1d(i)
    155         enddo
    156         close(11)
    157 
    158 ! extension pour klev+1 a nlaykim_tot avec tcp.ver
    159 ! la jonction klev/klev+1 est brutale... Tant pis ?
    160         ialt=1
    161         do l=klev+1,nlaykim_tot
    162            zalt=rmil(l)-rad/1000.
    163            do i=ialt,130
    164             if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
    165               ialt=i
    166             endif
    167            enddo
    168            factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
    169            press_c(l) = exp(  log(p1d(ialt))   * (1-factalt)  &
    170                            + log(p1d(ialt+1)) * factalt    )
    171            nb(l)      = exp(  log(ct1d(ialt))   * (1-factalt) & 
    172                            + log(ct1d(ialt+1)) * factalt    )
    173            temp_c(l)  = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
    174 !          print*,l,zalt,press_c(l),nb(l),temp_c(l)
    175         enddo
    176        
    177 ! caracteristiques des composes:       
    178 ! -----------------------------
    179         mass(:) = 0.0
    180         call comp(nomqy_c,nb,temp_c,mass,md)
    181         print*,'           Mass'
    182         do ic=1,nkim
    183           print*,nomqy_c(ic),mass(ic)
    184 !         if(nomqy_c(ic).eq.'CH4') print*,"MD=",md(:,ic)
    185         enddo   
    186                
    187 ! taux de photodissociations:
    188 ! --------------------------
    189         call disso(krpd,nbp_lat)
    190 
    191 ! reactions chimiques:
    192 ! -------------------
    193         call chimie(nomqy_c,nb,temp_c,krate,reactif, &
    194                     nom_perte,nom_prod,perte,prod)
    195 !        print*,'nom_prod, nom_perte:'
    196 !        do ic=1,nkim
    197 !          print*,nom_prod(ic),nom_perte(ic)
    198 !        enddo
    199 !        print*,'premieres prod, perte(1:reaction,2:compagnon):'
    200 !        do ic=1,nkim
    201 !          print*,prod(1,ic),perte(1,1,ic),perte(2,1,ic)
    202 !        enddo
    203 
    204 !       l = klev-3
    205 !       print*,'krate a p=',press_c(l),' reactifs et produits:'
    206 !       do ic=1,nd_kim+1
    207 !         print*,ic,krpd(7,ic,l,4)*nb(l),"  ",
    208 !    .     nomqy_c(reactif(1,ic)+1),
    209 !    .     nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
    210 !    .     nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
    211 !       enddo
    212 !       do ic=nd_kim+2,nr_kim
    213 !         print*,ic,krate(l,ic),"  ",
    214 !    .     nomqy_c(reactif(1,ic)+1),
    215 !    .     nomqy_c(reactif(2,ic)+1),nomqy_c(reactif(3,ic)+1),
    216 !    .     nomqy_c(reactif(4,ic)+1),nomqy_c(reactif(5,ic)+1)
    217 !       enddo
    218 
    219 
    220 !   init kedd
    221 !   ---------
    222 !   kedd stays constant with time and space
    223 !   => linked to pressure rather than altitude...
     364           rinter(l)  = (czlev(ig,l)+rad)/1000.0             ! km
     365           rmil(l)    = (czlay(ig,l)+rad)/1000.0             ! km
     366           temp_c(l)  = ctemp(ig,l)                          ! K
     367           phi_c(l)   = cpphi(ig,l)                          ! m2.s-2
     368           press_c(l) = cplay(ig,l)/100.                     ! mbar
     369           nb(l)      = 1.e-4*press_c(l) / (kbol*temp_c(l))  ! cm-3
     370        ENDDO
     371        rinter(klev+1)=(czlev(ig,klev+1)+rad)/1000.
     372
     373        ! b. Extension in upper atmosphere with Vervack profile
     374        !    NB : The transition klev/klev+1 is maybe harsh if we have correct g above but not under 
     375
     376        ipres=1
     377        DO l=klev+1,nlaykim_tot
     378          press_c(l) = preskim(l-klev) / 100.0
     379          DO i=ipres,130
     380              IF ( (press_c(l).LE.p1d(i)) .AND. (press_c(l).GT.p1d(i+1)) ) THEN
     381                ipres=i
     382              ENDIF
     383          ENDDO
     384          factp = (press_c(l)-p1d(ipres)) / (p1d(ipres+1)-p1d(ipres))
     385
     386          nb(l)      = exp( log(ct1d(ipres))*(1-factp) + log(ct1d(ipres+1))* factp )
     387          temp_c(l)  = t1d(ipres)*(1-factp) + t1d(ipres+1)*factp
     388        ENDDO
    224389 
    225       kedd(:) = 5e2 ! valeur mise par defaut 
    226                ! UNITE ? doit etre ok pour gptitan
    227       do l=1,nlaykim_tot
    228        zalt=rmil(l)-rad/1000.  ! z en km
    229        if     ((zalt.ge.250.).and.(zalt.le.400.)) then
    230          kedd(l) = 10.**(3.+(zalt-250.)/50.)
    231 ! 1E3 at 250 km
    232 ! 1E6 at 400 km
    233        elseif ((zalt.gt.400.).and.(zalt.le.600.)) then
    234          kedd(l) = 10.**(6.+1.3*(zalt-400.)/200.)
    235 ! 2E7 at 600 km
    236        elseif ((zalt.gt.600.).and.(zalt.le.900.)) then
    237          kedd(l) = 10.**(7.3+0.7*(zalt-600.)/300.)
    238 ! 1E8 above 900 km
    239        elseif ( zalt.gt.900.                    ) then
    240         kedd(l) = 1.e8
    241        endif
    242       enddo
    243 
    244       ENDIF  ! premier appel
    245 
    246 !***********************************************************************
    247 !-----------------------------------------------------------------------
    248 
    249 !   calcul declin_c (en degres)
    250 !   ---------------------------
    251 
    252       declin_c = declin_rad*180./pi
    253 !      print*,'declinaison en degre=',declin_c
    254        
    255 !***********************************************************************
    256 !***********************************************************************
    257 !
    258 !                BOUCLE SUR LES LATITUDES
    259 !
    260 !                * Permet de faire le calcul une seule fois par lat
    261 !
    262       DO j=1,ngrid
    263      
    264       if (j.eq.1) then
    265          jm1=1
    266       else
    267          jm1=j-1
    268       endif
    269 
    270       if((j.eq.1).or.(klat(j).ne.klat(jm1))) then
    271 
    272 !***********************************************************************
    273 !***********************************************************************
    274 
    275 !-----------------------------------------------------------------------
    276 !
    277 !   Distance radiale (intercouches, en km)
    278 !   ----------------------------------------
    279 
    280        do l=1,klev
    281          rinter(l) = (rad+czlev(j,l))/1000.
    282          rmil(l)   = (rad+czlay(j,l))/1000.
    283 !        print*,rinter(l)
    284        enddo
    285        rinter(klev+1)=(rad+czlev(j,klev+1))/1000.
    286 
    287 ! au-dessus du GCM, dr regulier et rinter(nlaykim_tot)=1290+2575 km.
    288        do l=klev+2,nlaykim_tot
    289          rinter(l) = rinter(klev+1)  &
    290                + (l-klev-1)*(3865.-rinter(klev+1))/(nlaykim_tot-klev-1)
    291          rmil(l-1) = (rinter(l-1)+rinter(l))/2.
    292        enddo
    293        rmil(nlaykim_tot) = rinter(nlaykim_tot)+(rinter(nlaykim_tot)-rinter(nlaykim_tot-1))/2.
    294 
    295 !-----------------------------------------------------------------------
    296 !
    297 !   Temperature, pression (mbar), densite (cm-3)
    298 !   -------------------------------------------
    299 
    300        DO l=1,klev
    301 !     temp_c (K):
    302                temp_c(l)  = ctemp(j,l)
    303 !     press_c (mbar):
    304                press_c(l) = cplay(j,l)/100.
    305 !     nb (cm-3):
    306                nb(l) = 1.e-4*press_c(l) / (kbol*temp_c(l))
    307        ENDDO
    308 ! extension pour klev+1 a nlaykim_tot avec tcp.ver
    309        ialt=1
    310        do l=klev+1,nlaykim_tot
    311            zalt=rmil(l)-rad/1000.
    312            do i=ialt,130
    313             if ((zalt.ge.r1d(i)).and.(zalt.lt.r1d(i+1))) then
    314               ialt=i
    315             endif
    316            enddo
    317            factalt = (zalt-r1d(ialt))/(r1d(ialt+1)-r1d(ialt))
    318            press_c(l) = exp(  log(p1d(ialt))   * (1-factalt) &
    319                            + log(p1d(ialt+1)) * factalt    )
    320            nb(l)      = exp(  log(ct1d(ialt))   * (1-factalt) &
    321                            + log(ct1d(ialt+1)) * factalt    )
    322            temp_c(l)  = t1d(ialt) * (1-factalt) + t1d(ialt+1) * factalt
    323        enddo
    324                
    325 !-----------------------------------------------------------------------
    326 !
    327 !   passage krpd => krate
    328 !   ---------------------
    329 ! initialisation krate pour dissociations
    330 
    331       if ((declin_c*10+267).lt.14.) then
    332           declinint = 0
    333           dec       = 0
    334       else
    335        if ((declin_c*10+267).gt.520.) then
    336           declinint = 14
    337           dec       = 534
    338        else
    339           declinint = 1
    340           dec       = 27
    341           do while( (declin_c*10+267).ge.real(dec+20) )
    342             dec       = dec+40
    343             declinint = declinint+1
    344           enddo
    345        endif
    346       endif
    347       if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
    348           factdec = ( declin_c - (dec-267)/10. ) / 4.
    349       else
    350           factdec = ( declin_c - (dec-267)/10. ) / 2.7
    351       endif
    352 
    353       do l=1,nlaykim_tot
    354 
    355 ! INTERPOL EN ALT POUR k (krpd tous les deux km)
    356         ialt    = int((rmil(l)-rad/1000.)/2.)+1
    357         factalt = (rmil(l)-rad/1000.)/2.-(ialt-1)
    358 
    359         do i=1,nd_kim+1
    360           krpddecm1 = krpd(declinint  ,i,ialt  ,klat(j)) * (1-factalt) &
    361                    + krpd(declinint  ,i,ialt+1,klat(j)) * factalt
    362           krpddec   = krpd(declinint+1,i,ialt  ,klat(j)) * (1-factalt) &
    363                    + krpd(declinint+1,i,ialt+1,klat(j)) * factalt
    364           krpddecp1 = krpd(declinint+2,i,ialt  ,klat(j)) * (1-factalt) &
    365                    + krpd(declinint+2,i,ialt+1,klat(j)) * factalt
    366 
    367                     ! nd_kim+1 correspond a la dissociation de N2 par les GCR
    368           if ( factdec.lt.0. ) then
    369         krate(l,i) = krpddecm1 * abs(factdec)  &
    370                   + krpddec   * ( 1 + factdec)
    371           endif
    372           if ( factdec.gt.0. ) then
    373         krate(l,i) = krpddecp1 * factdec       &
    374                   + krpddec   * ( 1 - factdec)
    375           endif
    376           if ( factdec.eq.0. ) krate(l,i) = krpddec
    377         enddo       
    378       enddo       
    379 
    380 !-----------------------------------------------------------------------
    381 !
    382 !   composition
    383 !   ------------
    384 
    385        do ic=1,nkim
    386         do l=1,klev
    387           cqy(l,ic)  = qy_c(j,l,ic)
    388           cqy0(l,ic) = cqy(l,ic)
    389         enddo
    390        enddo
    391 
    392 ! lecture du fichier compo_klat(j) (01 à 49) pour klev+1 a nlaykim_tot
    393 
    394       WRITE(str2,'(i2.2)') klat(j)
    395       fich = "comp_"//str2//".dat"
    396       inquire (file=fich,exist=okfic)
    397       if (okfic) then
    398        open(11,file=fich,status='old')
    399 ! premiere ligne=declin
    400        read(11,'(A15)') ficdec
    401        write(curdec,'(E15.5)') declin_c
    402 ! si la declin est la meme: ce fichier a deja ete reecrit
    403 ! on lit la colonne t-1 au lieu de la colonne t
    404 ! (cas d une bande de latitude partagee par 2 procs)
    405        do ic=1,nkim
    406         read(11,'(A10)') name
    407         if (name.ne.nomqy_c(ic)) then
    408           print*,"probleme lecture ",fich
    409           print*,name,nomqy_c(ic)
    410           stop
     390        ! We build altitude with hydrostatic equilibrium on preskim grid with Vervack profile
     391        ! ( keeping in mind that preskim is built based on Vervack profile with dz=10km )
     392
     393        DO l=klev+1,nlaykim_tot
     394
     395           ! Compute geopotential on the upper grid, to have correct altitude
     396           ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     397           ! NB : Upper chemistry model is based on observational profiles and needs correct altitudes !!
     398           ! Altitudes can (and must) be calculated with correct g even if it's not the case in the GCM under
     399
     400           temp1 = 0.5*(temp_c(l-1)+temp_c(l)) ! interlayer temp
     401           phi_c(l) = phi_c(l-1) + r*temp1*log(press_c(l-1)/press_c(l)) ! Geopotential assuming hydrostatic equilibrium
     402
     403           rmil(l)  =  ( g*rad*rad / (g*rad - phi_c(l)) ) / 1000.0 ! z(phi) with g varying with altitude
     404        ENDDO
     405
     406        DO l=klev+2,nlaykim_tot
     407          rinter(l) = 0.5*(rmil(l-1) + rmil(l))
     408        ENDDO
     409
     410        ! 2. From krpd, compute krate for dissociations (declination-latitude-altitude interpolation)
     411        ! -------------------------------------------------------------------------------------------
     412
     413        ! a. Calculate declination dependence
     414
     415        if ((declin_c*10+267).lt.14.) then
     416           idec = 0
     417           dec  = 0
     418        else
     419           if ((declin_c*10+267).gt.520.) then
     420              idec = 14
     421              dec  = 534
     422           else
     423              idec = 1
     424              dec  = 27
     425              do while( (declin_c*10+267).ge.real(dec+20) )
     426                 dec  = dec+40
     427                 idec = idec+1
     428              enddo
     429           endif
    411430        endif
    412         if (ficdec.eq.curdec) then
    413           do l=klev+1,nlaykim_tot
    414             read(11,*) ficalt,cqy(l,ic),newq
    415           enddo
     431        if ((declin_c.ge.-24.).and.(declin_c.le.24.)) then
     432           factdec = ( declin_c - (dec-267)/10. ) / 4.
    416433        else
    417           do l=klev+1,nlaykim_tot
    418             read(11,*) ficalt,oldq,cqy(l,ic)
    419           enddo
     434           factdec = ( declin_c - (dec-267)/10. ) / 2.7
    420435        endif
    421        enddo
    422        close(11)
    423       else       ! le fichier n'est pas la
    424        do ic=1,nkim
    425         do l=klev+1,nlaykim_tot
    426           cqy(l,ic)=cqy(klev,ic)
    427         enddo
    428        enddo
    429       endif
    430       cqy0 = cqy
    431 
    432 !-----------------------------------------------------------------------
    433 !
    434 !   Appel de chimietitan
    435 !   --------------------
    436        
    437        call gptitan(rinter,temp_c,nb,                          &
    438                    nomqy_c,cqy,                               &
    439                    duree,(klat(j)-1),mass,md,                 &
    440                    kedd,botCH4,krate,reactif,                 &
    441                    nom_prod,nom_perte,prod,perte,             &
    442                    aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,  &
    443                    htoh2,surfhaze)
    444        
    445 !   Tendances composition
    446 !   ---------------------
    447      
    448        do ic=1,nkim
    449         do l=1,klev
    450           dqyc(j,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim  ! en /s
    451         enddo
    452        enddo
    453 
    454 
    455 !-----------------------------------------------------------------------
    456 !
    457 !   sauvegarde compo sur nlaykim_tot
    458 !   -------------------------
    459 
    460 ! dans fichier compo_klat(j) (01 à 49)
    461        
    462       open(11,file=fich,status='replace')
    463 ! premiere ligne=declin
    464       write(11,'(E15.5)') declin_c
    465       do ic=1,nkim
    466        write(11,'(A10)') nomqy_c(ic)
    467        do l=klev+1,nlaykim_tot
    468         write(11,'(E15.5,1X,E15.5,1X,E15.5)') rmil(l)-rad/1000.,cqy0(l,ic),cqy(l,ic)
    469        enddo
    470       enddo
    471       close(11)
    472 
    473 !***********************************************************************
    474 !***********************************************************************
    475 !              FIN: BOUCLE SUR LES LATITUDES
    476 
    477       else      ! same latitude, we don't do calculations again, only adjust longitudinal variations
    478         dqyc(j,:,:) = dqyc(jm1,:,:)/qy_c(jm1,:,:)*qy_c(j,:,:)
    479       endif
    480 
    481       ENDDO
    482      
    483 !***********************************************************************
    484 !***********************************************************************
    485 
    486       firstcall = .false.
    487      
    488       end subroutine calchim
     436
     437        ! b. Calculate klat for interpolation on fixed latitudes of actinic fluxes input
     438
     439        klat=1
     440        DO i=1,nlat_actfluxes
     441          IF (latitude(ig).LT.lat_actfluxes(i)) klat=i
     442        ENDDO
     443        IF (klat==nlat_actfluxes) THEN ! avoid rounding problems
     444          klat    = nlat_actfluxes-1
     445          factlat = 1.0
     446        ELSE
     447          factlat = (latitude(ig)-lat_actfluxes(klat))/(lat_actfluxes(klat+1)-lat_actfluxes(klat))
     448        ENDIF
     449
     450        ! c. Altitude loop
     451
     452        DO l=1,nlaykim_tot
     453
     454           ! Calculate ialt for interpolation in altitude (krpd every 2 km)
     455           ialt    = int((rmil(l)-rad/1000.)/2.)+1
     456           factalt = (rmil(l)-rad/1000.)/2.-(ialt-1)
     457
     458           ! Altitude can go above top limit of UV levels - in this case we keep the 1310km top fluxes
     459           IF (ialt.GT.nlrt_kim-1) THEN
     460             ialt     = nlrt_kim-1 ! avoid out-of-bound array
     461             factalt  = 1.0
     462           ENDIF
     463
     464           DO i=1,nd_kim+1
     465
     466              krpddecm1 =   (   krpd(idec  ,i,ialt  ,klat)   * (1.0-factalt)                   &
     467                              + krpd(idec  ,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
     468                          * (   krpd(idec  ,i,ialt  ,klat+1) * (1.0-factalt)                   &
     469                              + krpd(idec  ,i,ialt+1,klat+1) * factalt       ) * factlat
     470              krpddec   =   (   krpd(idec+1,i,ialt  ,klat)   * (1.0-factalt)                   &
     471                              + krpd(idec+1,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
     472                          * (   krpd(idec+1,i,ialt  ,klat+1) * (1.0-factalt)                   &
     473                              + krpd(idec+1,i,ialt+1,klat+1) * factalt       ) * factlat
     474              krpddecp1 =   (   krpd(idec+2,i,ialt  ,klat)   * (1.0-factalt)                   &
     475                              + krpd(idec+2,i,ialt+1,klat)   * factalt       ) * (1.0-factlat) &
     476                          * (   krpd(idec+2,i,ialt  ,klat+1) * (1.0-factalt)                   &
     477                              + krpd(idec+2,i,ialt+1,klat+1) * factalt       ) * factlat
     478
     479              ! nd_kim+1 is dissociation of N2 by GCR
     480              if ( factdec.lt.0. ) then
     481                 krate(l,i) = krpddecm1 * abs(factdec) + krpddec   * ( 1.0 + factdec)
     482              endif
     483              if ( factdec.gt.0. ) then
     484                 krate(l,i) = krpddecp1 * factdec      + krpddec   * ( 1.0 - factdec)
     485              endif
     486              if ( factdec.eq.0. ) krate(l,i) = krpddec
     487
     488           ENDDO ! i=1,nd_kim+1
     489        ENDDO ! l=1,nlaykim_tot
     490
     491        ! 3. Read composition
     492        ! -------------------
     493
     494        DO ic=1,nkim
     495           DO l=1,klev
     496              cqy(l,ic) = qy_c(ig,l,ic) ! advected tracers for the GCM part converted to molar frac.
     497           ENDDO
     498           
     499           DO l=1,nlaykim_up
     500              cqy(klev+l,ic) = ykim_up(ic,ig,l) ! ykim_up for the upper atm.
     501           ENDDO
     502        ENDDO
     503
     504        cqy0(:,:) = cqy(:,:) ! Stores compo. before modifs
     505
     506        ! 4. Call main Titan chemistry C routine
     507        ! --------------------------------------
     508
     509        call gptitan(rinter,temp_c,nb,                  &
     510             nomqy_c,cqy,                               &
     511             dtchim,latitude(ig)*180./pi,mass,md,       &
     512             kedd,botCH4,krate,reactif,                 &
     513             nom_prod,nom_perte,prod,perte,             &
     514             aerprod,utilaer,cmaer,cprodaer,ccsn,ccsh,  &
     515             htoh2,surfhaze)
     516
     517        ! 5. Calculates tendencies on composition for advected tracers
     518        ! ------------------------------------------------------------
     519        DO ic=1,nkim
     520           DO l=1,klev
     521              dqyc(ig,l,ic) = (cqy(l,ic) - cqy0(l,ic))/dtchim ! (mol/mol/s)
     522           ENDDO
     523        ENDDO
     524
     525        ! 6. Update ykim_up
     526        ! -----------------
     527        DO ic=1,nkim
     528           DO l=1,nlaykim_up
     529              ykim_up(ic,ig,l) = cqy(klev+l,ic)
     530           ENDDO
     531        ENDDO
     532        ! NB: The full vertical composition grid will be created only for the outputs
     533
     534
     535     ELSE ! In 2D chemistry, if following grid point at same latitude, same zonal mean so don't do calculations again !
     536        dqyc(ig,:,:)    = dqyc(igm1,:,:) ! will be put back in 3D with longitudinal variations assuming same relative tendencies within a lat band
     537        ykim_up(:,ig,:) = ykim_up(:,igm1,:) ! no horizontal mixing in upper layers -> no longitudinal variations
     538     ENDIF
     539
     540  ENDDO
     541
     542END SUBROUTINE calchim

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

    r1904 r1908  
    11
    22
    3 SUBROUTINE calmufi(dt, plev, zlev, play, zlay, temp, pq, zdq)
     3SUBROUTINE calmufi(dt, plev, zlev, play, zlay, temp, pq, zdqfi, zdq)
    44  !! Interface subroutine to YAMMS model for Titan LMDZ GCM.
    55  !!
    6   !! The subroutine computes the microphysics processes for a si:le vertical column.
     6  !! The subroutine computes the microphysics processes for a single vertical column.
    77  !!
    88  !! - All input vectors are assumed to be defined from GROUND to TOP of the atmosphere.
    … …  
    3030  REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp  !! Temperature at the center of each layer (K).
    3131
    32   REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq   !! Tracers (\(kg.kg^{-1}}\)).
    33   REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq  !! Tendency for tracers (\(kg.kg^{-1}}\)).
     32  REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq    !! Tracers (\(kg.kg^{-1}}\)).
     33  REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqfi !! Tendency from former processes for tracers (\(kg.kg^{-1}}\)).
     34  REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq   !! Microphysical tendency for tracers (\(kg.kg^{-1}}\)).
     35 
     36  REAL(kind=8), DIMENSION(:,:,:), ALLOCATABLE :: zq !! Local tracers updated from former processes (\(kg.kg^{-1}}\)).
    3437 
    3538  REAL(kind=8), DIMENSION(:), ALLOCATABLE :: m0as !! 0th order moment of the spherical mode (\(m^{-2}\)).
    … …  
    5659
    5760  INTEGER :: ilon, i,nices
    58   INTEGER :: nlon,nlay
     61  INTEGER :: nq,nlon,nlay
    5962
    6063  ! Read size of arrays
     64  nq    = size(pq,DIM=3)
    6165  nlon  = size(play,DIM=1)
    6266  nlay  = size(play,DIM=2)
    … …  
    6569  ALLOCATE( int2ext(nlay) ) 
    6670
    67   ! Loop on horizontal grid points
     71  ! Allocate arrays
     72  ALLOCATE( zq(nlon,nlay,nq) )
    6873
    69   ! Allocate arrays
    7074  ALLOCATE( m0as(nlay) )
    7175  ALLOCATE( m3as(nlay) )
    … …  
    8993  zdq(:,:,:) = 0.0
    9094
     95  ! Initialize tracers updated with former processes from physics
     96  zq(:,:,:) = pq(:,:,:) + zdqfi(:,:,:)*dt
     97
     98  ! Loop on horizontal grid points
    9199  DO ilon = 1, nlon
    92100    ! Convert tracers to extensive ( except for gazs where we work with molar mass ratio )
    93101    ! We suppose a given order of tracers !
    94102    int2ext(:) = ( plev(ilon,1:nlay)-plev(ilon,2:nlay+1) ) / g
    95     m0as(:) = pq(ilon,:,1) * int2ext(:)
    96     m3as(:) = pq(ilon,:,2) * int2ext(:)
    97     m0af(:) = pq(ilon,:,3) * int2ext(:)
    98     m3af(:) = pq(ilon,:,4) * int2ext(:)
     103    m0as(:) = zq(ilon,:,1) * int2ext(:)
     104    m3as(:) = zq(ilon,:,2) * int2ext(:)
     105    m0af(:) = zq(ilon,:,3) * int2ext(:)
     106    m3af(:) = zq(ilon,:,4) * int2ext(:)
    99107   
    100108    if (callclouds) then ! if call clouds
    101       m0n(:) = pq(ilon,:,5) * int2ext(:)
    102       m3n(:) = pq(ilon,:,6) * int2ext(:)
     109      m0n(:) = zq(ilon,:,5) * int2ext(:)
     110      m3n(:) = zq(ilon,:,6) * int2ext(:)
    103111      do i=1,nices
    104         m3i(:,nices) = pq(ilon,:,6+i) * int2ext(:)
    105         gazs(:,i)    = pq(ilon,:,ices_indx(i)) / rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
     112        m3i(:,nices) = zq(ilon,:,6+i) * int2ext(:)
     113        gazs(:,i)    = zq(ilon,:,ices_indx(i)) / rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
    106114        ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one
    107115      enddo

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

    r1903 r1908  
    11MODULE comchem_h
    22
    3 ! ----------------------------------------------------------------------------
    4 ! Purpose : Stores data relative to chemistry in the GCM and upper chemistry
    5 ! -------
    6 !           Please note that upper fields ykim_up are in molar fraction !
     3! -------------------------------------------------------------------------------------------------------------
     4! Purpose : + Stores data relative to :  * 1. Chemistry in the GCM
     5! -------                                * 2. Upper chemistry pressure grid
     6!                                        * 3. Coupling with C photochem. module ( cf calchim.F90)
     7!           
     8!           + Also contains a routine of initialization for chemistry in the GCM.
    79!
    8 !          NB : For newstart there is a specific comchem_newstart_h module.
     10!           + NB : For newstart there is a specific comchem_newstart_h module.
    911!
    1012! Author : Jan Vatant d'Ollone (2017-18)
    … …  
    1214!
    1315! NB : A given order is assumed for the 44 chemistry tracers :
    14 ! --   H, H2, CH, CH2s, CH2, CH3, CH4, C2, C2H, C2H2, C2H3, C2H4, C2H5,
     16!      H, H2, CH, CH2s, CH2, CH3, CH4, C2, C2H, C2H2, C2H3, C2H4, C2H5,
    1517!      C2H6, C3H3, C3H5, C3H6, C3H7, C4H, C4H3, C4H4, C4H2s, CH2CCH2,
    1618!      CH3CCH, C3H8, C4H2, C4H6, C4H10, AC6H6, C3H2, C4H5, AC6H5, N2,
    1719!      N4S, CN, HCN, H2CN, CHCN, CH2CN, CH3CN, C3N, HC3N, NCCN, C4N2
    18        
    19 ! ----------------------------------------------------------------------------
     20!
     21!       
     22! IMPORTANT : Advected chem. tracers are in MASS fraction but upper fields ykim_up are in MOLAR fraction !
     23!
     24! -------------------------------------------------------------------------------------------------------------
    2025
    2126IMPLICIT NONE 
     27
     28   ! ~~~~~~~~~~~~~~~~~~~~~~~~
     29   ! 1. Chemistry in the GCM
     30   ! ~~~~~~~~~~~~~~~~~~~~~~~~
    2231
    2332   !! Hard-coded number of chemical species for Titan chemistry
    … …  
    5059       40.07  , 40.07 , 44.11, 50.06, 54.09, 58.13, 78.1136, 38.05, 53.07, 77.1136, 28.0134, &
    5160       14.01  , 26.02 , 27.04, 28.05, 39.05, 40.04, 41.05  , 50.04, 51.05, 52.04  , 76.1   /)
     61
     62   !! Hard-coded molar fraction of surface methane
     63   REAL, PARAMETER :: botCH4 = 0.048
    5264   
    53    ! !!!!!!!!!!!!!!!!!!!!
    54    !  Upper chemistry
    55    ! !!!!!!!!!!!!!!!!!!!!
    56    
    57    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    58    ! These parameters as well as nkim above, MUST match titan.h in chimtitan
    59    INTEGER, PARAMETER :: nd_kim   = 54     ! Number of photodissociations
    60    INTEGER, PARAMETER :: nr_kim   = 377    ! Number of reactions in chemistry scheme
    61    INTEGER, PARAMETER :: nlrt_kim = 650    ! For the UV rad. transf., 650 levels of 2 km
    62    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    63    
     65
     66
     67   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
     68   !  2. Upper chemistry grid
     69   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~
     70
    6471   INTEGER, SAVE :: nlaykim_up   ! Number of upper atm. layers for chemistry from GCM top to 4.5E-5 Pa (1300km)
    6572   INTEGER, SAVE :: nlaykim_tot  ! Number of total layers for chemistry from surface to 4.5E-5 Pa (1300km)
    6673!$OMP_THREADPRIVATE(nlaykim_up,nlay_kim_tot)
    6774
    68   ! NB : For the startfile we use nlaykim_up grid (upper atm) and for outputs we use nlaykim_tot grid (all layers)
     75   ! NB : For the startfile we use nlaykim_up grid (upper atm) and for outputs we use nlaykim_tot grid (all layers)
     76 
     77   REAL*8, PARAMETER :: grkim_dz = 10.0 ! Vertical discretization of the upper chemsitry grid (km)
    6978
    70    REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: preskim  ! Pressure (Pa) of upper chemistry (mid)-layers
    71 !$OMP_THREADPRIVATE(preskim)
     79   REAL,SAVE,ALLOCATABLE,DIMENSION(:)     :: preskim  ! Pressure (Pa) of upper chemistry (mid)-layers
     80   REAL,SAVE,ALLOCATABLE,DIMENSION(:,:)   :: zlay_kim ! Altitude (km) of all chemistry (mid)-layers
     81   REAL,SAVE,ALLOCATABLE,DIMENSION(:,:,:) :: ykim_up  ! Upper chemistry fields (mol/mol)
     82!$OMP_THREADPRIVATE(preskim,zlay_kim,ykim_up)
    7283
    73    REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: zlay_kim  ! Altitude (km) of all chemistry (mid)-layers
    74 !$OMP_THREADPRIVATE(zlay_kim)
    7584
    76    REAL,SAVE,ALLOCATABLE,DIMENSION(:,:,:) :: ykim_up ! Upper chemistry fields (mol/mol)
    77 !$OMP_THREADPRIVATE(ykim_up)
     85
     86   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     87   !  3. Interface with photochemical module (cf calchim.F90)
     88   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     89   
     90   ! These 3 parameters as well as nkim above, MUST match titan.h in chimtitan !!
     91   INTEGER, PARAMETER :: nd_kim   = 54     ! Number of photodissociations
     92   INTEGER, PARAMETER :: nr_kim   = 377    ! Number of reactions in chemistry scheme
     93   INTEGER, PARAMETER :: nlrt_kim = 650    ! For the UV rad. transf., 650 levels of 2 km
     94
     95   !! Hardcoded latitude discretisation for actinic fluxes - MUST be coherent with disso.c input files !!
     96   INTEGER, PARAMETER                         :: nlat_actfluxes = 49
     97   REAL, DIMENSION(nlat_actfluxes), PARAMETER :: lat_actfluxes  = (/ &
     98      90.0 ,  86.25,  82.5 ,  78.75,  75.0 ,  71.25,  67.5 ,  63.75,  60.0 ,  56.25,  52.5 ,  48.75, &
     99      45.0 ,  41.25,  37.5 ,  33.75,  30.0 ,  26.25,  22.5 ,  18.75,  15.0 ,  11.25,   7.5 ,   3.75, &
     100       0.0 ,  -3.75,  -7.5 , -11.25, -15.0 , -18.75, -22.5 , -26.25, -30.0 , -33.75, -37.5 , -41.25, &
     101     -45.0 , -48.75, -52.5 , -56.25, -60.0 , -63.75, -67.5 , -71.25, -75.0 , -78.75, -82.5 , -86.25, &
     102     -90.0 /)
     103   
    78104
    79105CONTAINS

  • trunk/LMDZ.TITAN/libf/phytitan/dyn1d/rcm1d.F

    r1844 r1908  
    1212     &                     dtemisice
    1313      use comdiurn_h, only: sinlat, coslat, sinlon, coslon
     14      use comchem_h, only: nkim, nlaykim_up, ykim_up
    1415      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
    1516      use phyredem, only: physdem0,physdem1
    … …  
    2324     &                            nday, iphysiq
    2425      use callkeys_mod, only: tracer,check_cpp_match,rings_shadow,
    25      &                        specOLR,pceil
     26     &                        specOLR,pceil,callchim
    2627      USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff, sig,
    2728     &                       presnivs,pseudoalt,scaleheight
    … …  
    670671         ENDDO
    671672         
    672 
    673 
    674673         DO ilayer=1,nlayer
    675674!     zlay(ilayer)=-300.E+0 *r*log(play(ilayer)/plev(1))
    … …  
    762761       
    763762      endif ! end if tracer
     763     
     764! Initialize upper atmosphere & chemistry
     765! ---------------------------------------
     766
     767      ! Calculate the # of upper chemistry layers
     768      call gr_kim_vervack
     769     
     770      write(*,*)
     771      write(*,*) " With the compiled vertical grid we found :"
     772      write(*,*) " Number of upper chemistry layers =", nlaykim_up
     773       
     774      if (callchim) then
     775           
     776        if (.NOT.allocated(ykim_up)) then
     777          allocate(ykim_up(nkim,1,nlaykim_up))
     778        endif
     779     
     780        ! Inquire if we have upper chemistry input profiles
     781       
     782        do iq=1,nkim
     783       
     784          file_prof = "profile_"//trim(noms(iq))//"_up"
     785
     786          INQUIRE(file=trim(file_prof),exist=findprof)
     787       
     788          IF (findprof) THEN
     789            ! Read profile for upper chemistry specie.
     790            PRINT *, "Input file found for upper chemistry specie ",
     791     &                trim(noms(iq))//"_up"
     792            OPEN(11,file=trim(file_prof),status='old',form='formatted')
     793            PRINT *, "---> reading upper profile from input file ..."
     794            DO ilayer=1,nlaykim_up
     795              READ (11,*) ykim_up(iq,1,ilayer)
     796            ENDDO
     797          ELSE
     798            PRINT *, "No input profile found for ",
     799     &                trim(noms(iq))//"_up"
     800            PRINT *, "---> we set it to zero - that's bold !!"
     801          ENDIF
     802
     803        enddo ! nkim loop
     804       
     805      endif ! of if callchim
     806     
    764807
    765808c  Write a "startfi" file

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

    r1903 r1908  
    159159       enddo
    160160
    161        print*, 'inicondens end'
    162        
    163161      end subroutine inicondens
Note: See TracChangeset for help on using the changeset viewer.