Ignore:
Timestamp:
Oct 17, 2023, 9:40:14 AM (16 months ago)
Author:
slebonnois
Message:

BdeBatz? : Cleans microphysics and makes few corrections for physics

Location:
trunk/LMDZ.TITAN/libf/phytitan
Files:
18 edited

Legend:

Unmodified
Added
Removed

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

    r3083 r3090  
    2121!   Lcx = Condensation latente heat [J.kg-1]
    2222!
    23 !   Author(s)
    24 !   ---------
     23!   Author
     24!   ------
    2525!   B. de Batz de Trenquelléon (07/2023)
    2626!
    … …  
    2828!   ---------------
    2929!   We suppose a given order of tracers !
     30!   Tracers come from microphysics (ices_indx) :
    3031!   CH4  --> 1
    3132!   C2H2 --> 2
    3233!   C2H6 --> 3
    33 !   HCN  --> 4 (TO DO)
     34!   HCN  --> 4
    3435!====================================================================   
    3536

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

    r3083 r3090  
    77!     Average the Rayleigh scattering in each band, weighting the
    88!     average by the blackbody function at temperature tstellar.
    9 !     Works for an arbitrary mix of gases (currently CO2, N2 and
    10 !     H2 are possible).
     9!     Works for an arbitrary mix of gases (currently N2, H2
     10!     and CH4 are possible).
    1111!     
    1212!     Authors
    … …  
    6464      typeknown=.false.
    6565      do igas=1,ngasmx
    66         if(gfrac(igas,nivref).lt.1.e-2)then
     66        if(gfrac(igas,nivref).lt.1.e-2)then
    6767            print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
    68             'as its mixing ratio is less than 0.01.'
     68               'as its mixing ratio is less than 0.01.'
    6969            ! ignore variable gas in Rayleigh calculation
    7070            ! ignore gases of mixing ratio < 0.01 in Rayleigh calculation
    … …  
    7878               ! uses n=1.000132 from Optics, Fourth Edition.
    7979               ! was out by a factor 2!
    80             elseif(igas.eq.igas_CH4)then
    81                 ! For CH4 we use the formalism of exo_k (developed by Jeremy Leconte)
    82                 ! the relation between exo_k formalism and LMDZ formalism is : (tauconsti*tauvari)=scalep*sigma_exok/(gravity*molecular_mass)
    83                 ! Values are taken from Caldas (2019), equation 12 & appendix D
    84                 ! 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 is a constant
    85                 ! 4/9=(2/3)**2 is approximately ((n+1)/(n**2+2))**2
    86                 ! 1E24 = (1E6)**4 because wavelength is in micron
    87                 ! 16.04*1.6726*1E-27 is CH4 molecular mass
    88                 tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  * scalep             
     80            elseif(igas.eq.igas_CH4)then
     81               ! For CH4 we use the formalism of exo_k (developed by Jeremy Leconte)
     82               ! the relation between exo_k formalism and LMDZ formalism is : (tauconsti*tauvari)=scalep*sigma_exok/(gravity*molecular_mass)
     83               ! Values are taken from Caldas (2019), equation 12 & appendix D
     84               ! 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 is a constant
     85               ! 4/9=(2/3)**2 is approximately ((n+1)/(n**2+2))**2
     86               ! 1E24 = (1E6)**4 because wavelength is in micron
     87               ! 16.04*1.6726*1E-27 is CH4 molecular mass
     88               tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  * scalep             
    8989            else
    9090               print*,'Error in calc_rayleigh: Gas species not recognised!'

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

    r3083 r3090  
    1212  !     J. Vatant d'Ollone
    1313  !       -> 2017 : Adapt to new physics, move to F90 and compute every grid point
    14   !       -> 2018 : Update saturation profiles from recent litterature ( Vuitton 2019 )
     14  !       -> 2018 : Update saturation profiles from recent litterature (Vuitton 2019)
     15  !     Modified
     16  !     --------
     17  !     B. de Batz de Trenquelléon
     18  !       -> 2022 : Update saturation profiles and correct the pressure unit
    1519  !     ==============================================================================
    1620
    … …  
    4044  ALLOCATE(x(nlon,nlay))
    4145
    42   !     By default, ysat=1, meaning we do not condense
     46  ! By default, ysat = 1, meaning we do not condense
    4347  ysat(:,:,:) = 1.0
    4448
    45   !     Main loop
    46 
    47   do ic=1,nkim
    48 
    49      if(trim(cnames(ic)).eq."CH4") then
    50 
    51         !where (temp(:,:).lt.90.65)                                     
    52         !   ysat(:,:,ic) = 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / temp(:,:)                  &
    53         !        -  115352.19e0 ) / temp(:,:) + 4055.6016e0 ) / temp(:,:)                    &
    54         !        + 453.92414e0 ) / temp(:,:) ) / press(:,:) * 1013.25e0
    55         !elsewhere
    56         !   ysat(:,:,ic) = 10.0**(3.901408e0 - ( ( 154567.02e0 / temp(:,:) - 1598.8512e0 )   &
    57         !        / temp(:,:) + 437.54809e0 ) / temp(:,:)) / press(:,:) * 1013.25e0
    58         !endwhere
    59 
    60         ! Fray and Schmidt (2009)
    61         ysat(:,:,ic) = (1000. / press(:,:)) * exp(1.051e1 - 1.110e3/temp(:,:)              &
    62                        - 4.341e3/temp(:,:)**2 + 1.035e5/temp(:,:)**3 - 7.910e5/temp(:,:)**4)
    63 
    64         ! Forcing CH4 to 1.4% minimum               
    65         where (ysat(:,:,ic).lt.0.014) ysat(:,:,ic) = 0.014
    66 
    67      else if(trim(cnames(ic)).eq."C2H2") then
    68 
    69 !        ysat(:,:,ic) = 10.0**(6.09748e0-1644.1e0/temp(:,:)+7.42346e0                        &
    70 !             * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
    71 
    72          ysat(:,:,ic) = (1000. / press(:,:)) * exp(1.340e1 - 2.536e3/temp(:,:)) ! Fray and Schmidt (2009)
    73 
    74      else if(trim(cnames(ic)).eq."C2H4") then
    75 
    76         where (temp(:,:).lt.89.0) ! not far from Fray and Schmidt, 2009
    77            ysat(:,:,ic) = 10.0**(1.5477e0 + (1.0e0/temp(:,:) - 0.011e0)                     &
    78                 * (16537.0e0*(1.0e0/temp(:,:) - 0.011e0) - 1038.1e0))                       &
    79                 / press(:,:) * 1.01325e0 / 760.0
    80         elsewhere (temp(:,:).lt.104.0) ! not far from Fray and Schmidt, 2009
    81            ysat(:,:,ic) = 10.0**(8.724e0 - 901.6e0/(temp(:,:) - 2.555e0) )                  &
    82                 / press(:,:) * 1013.25e0 / 760.0
    83         elsewhere (temp(:,:).lt.120.0)
    84            ysat(:,:,ic) = 10.0**(50.79e0 - 1703.0e0/temp(:,:) - 17.141e0                    &
    85                 * alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0
    86         elsewhere (temp(:,:).lt.155.0)
    87            ysat(:,:,ic) = 10.0**(6.74756e0 - 585.0e0/(temp(:,:) - 18.16e0) )                &
    88                 / press(:,:) * 1013.25e0 / 760.0
    89         endwhere
    90 
    91      else if(trim(cnames(ic)).eq."C2H6") then
    92 
    93 !        where (temp(:,:).lt.90.)
    94 !           ysat(:,:,ic) = 10.0**(10.01e0-1085.0e0/(temp(:,:)-0.561e0) )                    &
    95 !                / press(:,:) * 1013.25e0 / 760.0e0
    96 !           ysat(:,:,ic) = 1.E5 * exp ( 15.11 - 2207./temp(:,:) - 24110./(temp(:,:)**2)     &
    97 !                        + 7.744E5/(temp(:,:)**3) - 1.161E7/(temp(:,:)**4)                  &
    98 !                        + 6.763E7/(temp(:,:)**5) ) / press(:,:) ! Fray and Schmidt (2009)
    99 !        elsewhere
    100 !           ysat(:,:,ic) = 10.0**(5.9366e0 - 1086.17e0/temp(:,:) + 3.83464e0                &
    101 !                * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
    102 !        endwhere
    103          ysat(:,:,ic) = (1000. / press(:,:)) * exp(1.511e1 - 2.207e3/temp(:,:)            &
    104                       - 2.411e4/temp(:,:)**2 + 7.744e5/temp(:,:)**3 - 1.161e7/temp(:,:)**4  &
    105                       + 6.763e7/temp(:,:)**5)
    106 
    107      else if((trim(cnames(ic)).eq."CH3CCH") .or. (trim(cnames(ic)).eq."CH2CCH2")) then
    108 
    109         ysat(:,:,ic) = 10.0**(2.8808e0 - 4.5e0*(249.9e0 - temp(:,:))                        &
     49  ! Main loop
     50  do ic = 1, nkim
     51     
     52      ! CH4 :
     53      !------
     54      if(trim(cnames(ic)).eq."CH4") then
     55         !where (temp(:,:).lt.90.65)                                     
     56         !   ysat(:,:,ic) = 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / temp(:,:)                 &
     57         !        -  115352.19e0 ) / temp(:,:) + 4055.6016e0 ) / temp(:,:)                   &
     58         !        + 453.92414e0 ) / temp(:,:) ) / press(:,:) * 1013.25e0
     59         !elsewhere
     60         !   ysat(:,:,ic) = 10.0**(3.901408e0 - ( (154567.02e0 / temp(:,:) - 1598.8512e0)    &
     61         !        / temp(:,:) + 437.54809e0 ) / temp(:,:)) / press(:,:) * 1013.25e0
     62         !endwhere
     63
     64         ! Fray and Schmidt (2009)
     65         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp(1.051e1 - 1.110e3/temp(:,:)                &
     66                        - 4.341e3/temp(:,:)**2 + 1.035e5/temp(:,:)**3 - 7.910e5/temp(:,:)**4)
     67
     68         ! Forcing CH4 to 1.4% minimum               
     69         where (ysat(:,:,ic).lt.0.014) ysat(:,:,ic) = 0.014
     70
     71      ! C2H2 :
     72      !-------
     73      else if(trim(cnames(ic)).eq."C2H2") then
     74         !ysat(:,:,ic) = 10.0**(6.09748e0-1644.1e0/temp(:,:)+7.42346e0                       &
     75         !     * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
     76
     77         ! Fray and Schmidt (2009)
     78         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp(1.340e1 - 2.536e3/temp(:,:))
     79
     80      ! C2H4 :
     81      !-------
     82      else if(trim(cnames(ic)).eq."C2H4") then
     83         ! not far from Fray and Schmidt (2009)
     84         where (temp(:,:).lt.89.0)
     85            ysat(:,:,ic) = 10.0**(1.5477e0 + (1.0e0/temp(:,:) - 0.011e0)                     &
     86                  * (16537.0e0*(1.0e0/temp(:,:) - 0.011e0) - 1038.1e0))                      &
     87                  / press(:,:) * 1.01325e0 / 760.0
     88         elsewhere (temp(:,:).lt.104.0)
     89            ysat(:,:,ic) = 10.0**(8.724e0 - 901.6e0/(temp(:,:) - 2.555e0) )                  &
     90                  / press(:,:) * 1013.25e0 / 760.0
     91         elsewhere (temp(:,:).lt.120.0)
     92            ysat(:,:,ic) = 10.0**(50.79e0 - 1703.0e0/temp(:,:) - 17.141e0                    &
     93                  * alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0
     94         elsewhere (temp(:,:).lt.155.0)
     95            ysat(:,:,ic) = 10.0**(6.74756e0 - 585.0e0/(temp(:,:) - 18.16e0) )                &
     96                  / press(:,:) * 1013.25e0 / 760.0
     97         endwhere
     98     
     99      ! C2H6 :
     100      !-------
     101      else if(trim(cnames(ic)).eq."C2H6") then
     102         where (temp(:,:).lt.90.)
     103            !ysat(:,:,ic) = 10.0**(10.01e0-1085.0e0/(temp(:,:)-0.561e0) )                    &
     104            !      / press(:,:) * 1013.25e0 / 760.0e0
     105            ! Fray and Schmidt (2009)
     106            ysat(:,:,ic) = (1.e3 / press(:,:)) * exp ( 15.11 - 2207./temp(:,:)               &
     107                           - 24110./(temp(:,:)**2) + 7.744E5/(temp(:,:)**3)                  &
     108                           - 1.161E7/(temp(:,:)**4) + 6.763E7/(temp(:,:)**5))
     109         elsewhere
     110            ysat(:,:,ic) = 10.0**(5.9366e0 - 1086.17e0/temp(:,:) + 3.83464e0                 &
     111                  * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
     112         endwhere
     113     
     114      ! CH3CCH :
     115      !---------
     116      else if((trim(cnames(ic)).eq."CH3CCH") .or. (trim(cnames(ic)).eq."CH2CCH2")) then
     117         ysat(:,:,ic) = 10.0**(2.8808e0 - 4.5e0*(249.9e0 - temp(:,:))                        &
    110118             / (1.15e0*temp(:,:) - 37.485e0) ) / press(:,:) * 1013.25e0 / 760.0e0
    111119
    112      else if(trim(cnames(ic)).eq."C3H6")  then
    113 
    114         ysat(:,:,ic) = 10.0**(7.4463e0 - 1028.5654e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
    115 
     120      ! C3H6 :
     121      !-------
     122      else if(trim(cnames(ic)).eq."C3H6")  then
     123         ysat(:,:,ic) = 10.0**(7.4463e0 - 1028.5654e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
     124
     125      ! C3H8 :
     126      !-------
    116127     else if(trim(cnames(ic)).eq."C3H8")  then
    117 
    118128        ysat(:,:,ic) = 10.0**(7.217e0 - 994.30251e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
    119129
    120      else if((trim(cnames(ic)).eq."C4H2").or.(trim(cnames(ic)).eq."C4H2s")) then
    121 
    122         ysat(:,:,ic) = 10.0**(96.26781e0 - 4651.872e0/temp(:,:) - 31.68595e0                &
     130      ! C4H2 :
     131      !-------
     132      else if((trim(cnames(ic)).eq."C4H2").or.(trim(cnames(ic)).eq."C4H2s")) then
     133         ysat(:,:,ic) = 10.0**(96.26781e0 - 4651.872e0/temp(:,:) - 31.68595e0                &
    123134             * alog10(temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
    124 
    125      else if(trim(cnames(ic)).eq."C4H4")  then
    126 
    127         ysat(:,:,ic)= 1.0e3 * exp(9.3898e0 - 2203.57/(temp(:,:)-43.15e0) ) / press(:,:)
    128 
    129      else if(trim(cnames(ic)).eq."C4H6")  then
    130 
    131         ysat(:,:,ic)= 10.0**(2.8808e0 - 4.6e0*(262.3e0 - temp(:,:))                         &
     135     
     136      ! C4H4 :
     137      !-------
     138      else if(trim(cnames(ic)).eq."C4H4")  then
     139         ysat(:,:,ic)= 1.0e3 * exp(9.3898e0 - 2203.57/(temp(:,:)-43.15e0) ) / press(:,:)
     140
     141      ! C4H6 :
     142      !-------
     143      else if(trim(cnames(ic)).eq."C4H6")  then
     144         ysat(:,:,ic)= 10.0**(2.8808e0 - 4.6e0*(262.3e0 - temp(:,:))                         &
    132145             /(1.15e0*temp(:,:) - 39.345e0) ) / press(:,:) * 1013.25e0 / 760.0e0
    133146
    134      else if(trim(cnames(ic)).eq."C4H10")  then
    135 
    136         ysat(:,:,ic)= 10.0**(8.446e0 - 1461.2e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
    137 
    138      else if(trim(cnames(ic)).eq."C6H2")  then
    139 
    140         ysat(:,:,ic)= 10.0**(4.666e0 - 4956e0/temp(:,:) + 25.845e0 *                        &
     147      ! C4H10 :
     148      !--------
     149      else if(trim(cnames(ic)).eq."C4H10")  then
     150         ysat(:,:,ic)= 10.0**(8.446e0 - 1461.2e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
     151
     152      ! C6H2 :
     153      !-------
     154      else if(trim(cnames(ic)).eq."C6H2")  then
     155         ysat(:,:,ic)= 10.0**(4.666e0 - 4956e0/temp(:,:) + 25.845e0 *                        &
    141156             alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
    142157
    143      else if(trim(cnames(ic)).eq."C8H2")  then
    144 
    145         ysat(:,:,ic)= 10.0**(3.95e0 - 6613e0/temp(:,:) + 35.055e0 *                         &
     158      ! C8H2 :
     159      !-------
     160      else if(trim(cnames(ic)).eq."C8H2")  then
     161         ysat(:,:,ic)= 10.0**(3.95e0 - 6613e0/temp(:,:) + 35.055e0 *                         &
    146162             alog10(1.0e3/temp(:,:)) ) / press(:,:) * 1013.25e0 / 760.0e0
    147163
    148      else if(trim(cnames(ic)).eq."AC6H6")  then
    149 
    150         !x = 1.0e0 - temp(:,:) / 562.2e0
    151         !ysat(:,:,ic)= 48.9e3 * exp( ( 1.33213 * x**1.5 - 6.98273 * x                        &
    152         !     - x**3 * (2.62863 + 3.33399 * x**3) ) * 562.2e0/temp(:,:) ) / press(:,:)
    153 
    154         ysat(:,:,ic) = 1.E5 * exp (17.35-5663./temp(:,:)) / press(:,:) ! Fray and Schmidt (2009)
    155        
    156      else if(trim(cnames(ic)).eq."HCN")  then
    157 
    158         !ysat(:,:,ic)= 10.0**(8.6165e0 - 1516.5e0/(temp(:,:) - 26.2e0) ) / press(:,:) * 1013.25e0 / 760.0e0
    159 
    160         ysat(:,:,ic) = (1000. / press(:,:)) * exp (1.393e1 - 3.624e3/temp(:,:)             &
    161                      - 1.325e5/temp(:,:)**2 + 6.314e6/temp(:,:)**3 - 1.128e8/temp(:,:)**4) ! Fray and Schmidt (2009)
    162 
    163      else if(trim(cnames(ic)).eq."CH3CN")  then
    164 
    165         ysat(:,:,ic)= 10.0**(8.458e0 - 1911.7e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
    166 
    167      else if(trim(cnames(ic)).eq."C2H3CN")  then
    168 
    169         ysat(:,:,ic)= 10.0**(9.3051e0 - 2782.21/(temp(:,:) - 51.15e0) ) / press(:,:) * 1013.25e0 / 760.0e0
    170 
    171      else if(trim(cnames(ic)).eq."NCCN")  then
    172 
    173         ysat(:,:,ic)=  10.0**(7.454e0 - 1832e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
    174 
    175      else if(trim(cnames(ic)).eq."HC3N")  then
    176 
    177         ysat(:,:,ic)= 10.0**(7.7446e0 - 1453.5609e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
    178 
    179      else if(trim(cnames(ic)).eq."C4N2")  then
    180 
    181         ysat(:,:,ic)= 10.0**(8.269e0 - 2155.0e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
    182 
    183      endif
     164      ! C6H6 :
     165      !-------
     166      else if(trim(cnames(ic)).eq."AC6H6")  then
     167         !x = 1.0e0 - temp(:,:) / 562.2e0
     168         !ysat(:,:,ic)= 48.9e3 * exp( ( 1.33213 * x**1.5 - 6.98273 * x                        &
     169         !     - x**3 * (2.62863 + 3.33399 * x**3) ) * 562.2e0/temp(:,:) ) / press(:,:)
     170
     171         ! Fray and Schmidt (2009)
     172         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp (17.35-5663./temp(:,:))
     173
     174      ! HCN :
     175      !-------
     176      else if(trim(cnames(ic)).eq."HCN")  then
     177         !ysat(:,:,ic)= 10.0**(8.6165e0 - 1516.5e0/(temp(:,:) - 26.2e0) ) / press(:,:) * 1013.25e0 / 760.0e0
     178         
     179         ! Fray and Schmidt (2009)
     180         ysat(:,:,ic) = (1.e3 / press(:,:)) * exp (1.393e1 - 3.624e3/temp(:,:)               &
     181                        - 1.325e5/temp(:,:)**2 + 6.314e6/temp(:,:)**3 - 1.128e8/temp(:,:)**4)
     182
     183      ! CH3CN :
     184      !--------
     185      else if(trim(cnames(ic)).eq."CH3CN")  then
     186         ysat(:,:,ic)= 10.0**(8.458e0 - 1911.7e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
     187   
     188      ! C2H3CN :
     189      !-------
     190      else if(trim(cnames(ic)).eq."C2H3CN")  then
     191         ysat(:,:,ic)= 10.0**(9.3051e0 - 2782.21/(temp(:,:) - 51.15e0) ) / press(:,:) * 1013.25e0 / 760.0e0
     192
     193      ! NCCN :
     194      !-------
     195      else if(trim(cnames(ic)).eq."NCCN")  then
     196         ysat(:,:,ic)=  10.0**(7.454e0 - 1832e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
     197   
     198      ! HC3N :
     199      !-------
     200      else if(trim(cnames(ic)).eq."HC3N")  then
     201         ysat(:,:,ic)= 10.0**(7.7446e0 - 1453.5609e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
     202   
     203      ! C4N2 :
     204      !-------
     205      else if(trim(cnames(ic)).eq."C4N2")  then
     206         ysat(:,:,ic)= 10.0**(8.269e0 - 2155.0e0/temp(:,:) ) / press(:,:) * 1013.25e0 / 760.0e0
     207
     208      endif
    184209  enddo
    185210

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

    r3083 r3090  
    3535!     Emmanuel 01/2001, Forget 09/2001
    3636!     Robin Wordsworth (2009)
    37 !     Jan Vatant d'Ollone (2018) -> corrk recombining case
    38 !     Bruno de Batz de Trenquelléon (2023) --> Optics of Titan's haze and coulds
     37!     
     38!     Modified
     39!     --------
     40!     Jan Vatant d'Ollone (2018)
     41!        --> corrk recombining case
     42!     B. de Batz de Trenquelléon (2023)
     43!        --> Titan's haze and could optics
     44!        --> clear/dark columns method
     45!        --> optical diagnostics
    3946!
    4047!==================================================================
    … …  
    8087      REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! IR optical thickness of layers within narrowbands for diags ().
    8188      REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! VI optical thickness of layers within narrowbands for diags ().
     89      ! Diagnostics
     90      REAL,INTENT(OUT) :: popthi(ngrid,nlayer,L_NSPECTI,5)  ! IR optical properties of haze within narrowbands (dtau,tau,k,w,g).
    8291      REAL,INTENT(OUT) :: popthv(ngrid,nlayer,L_NSPECTV,5)  ! VI optical properties of haze within narrowbands (dtau,tau,k,w,g).
     92      REAL,INTENT(OUT) :: poptci(ngrid,nlayer,L_NSPECTI,5)  ! IR optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g).
    8393      REAL,INTENT(OUT) :: poptcv(ngrid,nlayer,L_NSPECTV,5)  ! VI optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g).
    84       REAL,INTENT(OUT) :: popthi(ngrid,nlayer,L_NSPECTI,5)  ! IR optical properties of haze within narrowbands (dtau,tau,k,w,g).
    85       REAL,INTENT(OUT) :: poptci(ngrid,nlayer,L_NSPECTI,5)  ! IR optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g).
    8694     
    8795     
    … …  
    97105      ! Optical values for the optci/cv subroutines
    98106      REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
    99      
     107      ! IR
    100108      REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     109      REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
     110      REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     111      REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     112      ! VI
     113      REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     114      REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
     115      REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
     116      REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     117      REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     118     
     119      ! Temporary optical values for the optci/cv subroutines (clear column)
    101120      REAL*8 dtaui_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     121      REAL*8 dtauv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     122      REAL*8 tauv_cc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
     123      REAL*8 taucumi_cc(L_LEVELS,L_NSPECTI,L_NGAUSS)
     124      REAL*8 taucumv_cc(L_LEVELS,L_NSPECTV,L_NGAUSS)
     125      REAL*8 wbari_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     126      REAL*8 wbarv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     127      REAL*8 cosbi_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     128      REAL*8 cosbv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     129      ! Temporary optical values for the optci/cv subroutines (dark column)
    102130      REAL*8 dtaui_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    103       REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    104       REAL*8 dtauv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    105131      REAL*8 dtauv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    106      
    107       REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    108       REAL*8 cosbv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     132      REAL*8 tauv_dc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
     133      REAL*8 taucumi_dc(L_LEVELS,L_NSPECTI,L_NGAUSS)
     134      REAL*8 taucumv_dc(L_LEVELS,L_NSPECTV,L_NGAUSS)
     135      REAL*8 wbari_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     136      REAL*8 wbarv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     137      REAL*8 cosbi_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    109138      REAL*8 cosbv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    110       REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    111       REAL*8 cosbi_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    112       REAL*8 cosbi_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    113      
    114       REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    115       REAL*8 wbari_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    116       REAL*8 wbari_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    117       REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    118       REAL*8 wbarv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    119       REAL*8 wbarv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    120      
    121       REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
    122       REAL*8 tauv_cc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
    123       REAL*8 tauv_dc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
    124      
    125       REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
    126       REAL*8 taucumv_cc(L_LEVELS,L_NSPECTV,L_NGAUSS)
    127       REAL*8 taucumv_dc(L_LEVELS,L_NSPECTV,L_NGAUSS)
    128       REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
    129       REAL*8 taucumi_cc(L_LEVELS,L_NSPECTI,L_NGAUSS)
    130       REAL*8 taucumi_dc(L_LEVELS,L_NSPECTI,L_NGAUSS)
    131      
     139 
     140      ! Optical diagnostics
     141      REAL*8 popt_hazei(L_LEVELS,L_NSPECTI,3)
    132142      REAL*8 popt_hazev(L_LEVELS,L_NSPECTV,3)
     143      REAL*8 popt_cloudsi(L_LEVELS,L_NSPECTI,3)
     144      REAL*8 popt_cloudsv(L_LEVELS,L_NSPECTV,3)
     145      ! Temporary optical diagnostics (clear column)
     146      REAL*8 popt_hazei_cc(L_LEVELS,L_NSPECTI,3)
    133147      REAL*8 popt_hazev_cc(L_LEVELS,L_NSPECTV,3)
     148      REAL*8 popt_cloudsi_cc(L_LEVELS,L_NSPECTI,3)
     149      REAL*8 popt_cloudsv_cc(L_LEVELS,L_NSPECTV,3)
     150      ! Temporary optical diagnostics (dark column)
     151      REAL*8 popt_hazei_dc(L_LEVELS,L_NSPECTI,3)
    134152      REAL*8 popt_hazev_dc(L_LEVELS,L_NSPECTV,3)
    135       REAL*8 popt_hazei(L_LEVELS,L_NSPECTI,3)
    136       REAL*8 popt_hazei_cc(L_LEVELS,L_NSPECTI,3)
    137       REAL*8 popt_hazei_dc(L_LEVELS,L_NSPECTI,3)
    138 
    139       REAL*8 popt_cloudsv(L_LEVELS,L_NSPECTV,3)
    140       REAL*8 popt_cloudsv_cc(L_LEVELS,L_NSPECTV,3)
     153      REAL*8 popt_cloudsi_dc(L_LEVELS,L_NSPECTI,3)
    141154      REAL*8 popt_cloudsv_dc(L_LEVELS,L_NSPECTV,3)
    142       REAL*8 popt_cloudsi(L_LEVELS,L_NSPECTI,3)
    143       REAL*8 popt_cloudsi_cc(L_LEVELS,L_NSPECTI,3)
    144       REAL*8 popt_cloudsi_dc(L_LEVELS,L_NSPECTI,3)
     155
    145156
    146157      REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
    … …  
    546557         endwhere
    547558         ! Diagnostics for clouds :
    548          where (popt_cloudsv_cc(:,:,1) .lt. 1.d-30)
    549             popt_cloudsv_cc(:,:,1) = 1.d-30
    550          endwhere
    551          where (popt_cloudsv_dc(:,:,1) .lt. 1.d-30)
    552             popt_cloudsv_dc(:,:,1) = 1.d-30
    553          endwhere
    554          popt_cloudsv(:,:,1) = -log( (1-Fcloudy)*exp(-popt_cloudsv_cc(:,:,1)) + Fcloudy*exp(-popt_cloudsv_dc(:,:,1)) )
    555          popt_cloudsv(:,:,2) = ((popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2)) &
    556                              / (popt_cloudsv_cc(:,:,1) + popt_cloudsv_dc(:,:,1)))
    557          where ((popt_cloudsv_cc(:,:,2).gt.0.0D0) .or. (popt_cloudsv_dc(:,:,2).gt.0.0D0))
    558             popt_cloudsv(:,:,3) = (popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2)*popt_cloudsv_cc(:,:,3) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2)*popt_cloudsv_dc(:,:,3)) &
    559                              / (popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2))
    560          elsewhere
    561             popt_cloudsv(:,:,3) = 0.0D0
    562          endwhere
     559         if (callclouds) then
     560            where (popt_cloudsv_cc(:,:,1) .lt. 1.d-30)
     561               popt_cloudsv_cc(:,:,1) = 1.d-30
     562            endwhere
     563            where (popt_cloudsv_dc(:,:,1) .lt. 1.d-30)
     564               popt_cloudsv_dc(:,:,1) = 1.d-30
     565            endwhere
     566            popt_cloudsv(:,:,1) = -log( (1-Fcloudy)*exp(-popt_cloudsv_cc(:,:,1)) + Fcloudy*exp(-popt_cloudsv_dc(:,:,1)) )
     567            popt_cloudsv(:,:,2) = ((popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2)) &
     568                              / (popt_cloudsv_cc(:,:,1) + popt_cloudsv_dc(:,:,1)))
     569            where ((popt_cloudsv_cc(:,:,2).gt.0.0D0) .or. (popt_cloudsv_dc(:,:,2).gt.0.0D0))
     570               popt_cloudsv(:,:,3) = (popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2)*popt_cloudsv_cc(:,:,3) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2)*popt_cloudsv_dc(:,:,3)) &
     571                              / (popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2))
     572            elsewhere
     573               popt_cloudsv(:,:,3) = 0.0D0
     574            endwhere
     575         endif
    563576         
    564577         if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
    … …  
    666679         endwhere
    667680         ! Diagnostics for clouds :
    668          where (popt_cloudsi_cc(:,:,1) .lt. 1.d-30)
    669             popt_cloudsi_cc(:,:,1) = 1.d-30
    670          endwhere
    671          where (popt_cloudsi_dc(:,:,1) .lt. 1.d-30)
    672             popt_cloudsi_dc(:,:,1) = 1.d-30
    673          endwhere
    674          popt_cloudsi(:,:,1) = -log( (1-Fcloudy)*exp(-popt_cloudsi_cc(:,:,1)) + Fcloudy*exp(-popt_cloudsi_dc(:,:,1)) )
    675          popt_cloudsi(:,:,2) = ((popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2)) &
    676                              / (popt_cloudsi_cc(:,:,1) + popt_cloudsi_dc(:,:,1)))
    677          where ((popt_cloudsi_cc(:,:,2).gt.0.0D0) .or. (popt_cloudsi_dc(:,:,2).gt.0.0D0))
    678             popt_cloudsi(:,:,3) = (popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2)*popt_cloudsi_cc(:,:,3) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2)*popt_cloudsi_dc(:,:,3)) &
    679                              / (popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2))
    680          elsewhere
    681             popt_cloudsi(:,:,3) = 0.0D0
    682          endwhere
     681         if (callclouds) then
     682            where (popt_cloudsi_cc(:,:,1) .lt. 1.d-30)
     683               popt_cloudsi_cc(:,:,1) = 1.d-30
     684            endwhere
     685            where (popt_cloudsi_dc(:,:,1) .lt. 1.d-30)
     686               popt_cloudsi_dc(:,:,1) = 1.d-30
     687            endwhere
     688            popt_cloudsi(:,:,1) = -log( (1-Fcloudy)*exp(-popt_cloudsi_cc(:,:,1)) + Fcloudy*exp(-popt_cloudsi_dc(:,:,1)) )
     689            popt_cloudsi(:,:,2) = ((popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2)) &
     690                              / (popt_cloudsi_cc(:,:,1) + popt_cloudsi_dc(:,:,1)))
     691            where ((popt_cloudsi_cc(:,:,2).gt.0.0D0) .or. (popt_cloudsi_dc(:,:,2).gt.0.0D0))
     692               popt_cloudsi(:,:,3) = (popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2)*popt_cloudsi_cc(:,:,3) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2)*popt_cloudsi_dc(:,:,3)) &
     693                              / (popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2))
     694            elsewhere
     695               popt_cloudsi(:,:,3) = 0.0D0
     696            endwhere
     697         endif
    683698
    684699         call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
    … …  
    758773     
    759774     
    760       ! Optical properties of haze and clouds (dtau, tau, k, w, g) :
     775      ! Diagnostics : optical properties of haze and clouds (dtau, tau, k, w, g) :
    761776      do l = 2, L_LEVELS, 2
    762777         lev2lay = L_NLAYRAD+1 - l/2
    763          
    764778         ! Visible :
    765779         do nw = 1, L_NSPECTV
    … …  
    768782            ! Optical thickness (dtau) :
    769783            popthv(ig,lev2lay,nw,1) = (popt_hazev(l,nw,1) + popt_hazev(l+1,nw,1)) / 2.0
    770             poptcv(ig,lev2lay,nw,1) = (popt_cloudsv(l,nw,1) + popt_cloudsv(l+1,nw,1)) / 2.0
     784            if (callclouds) then
     785               poptcv(ig,lev2lay,nw,1) = (popt_cloudsv(l,nw,1) + popt_cloudsv(l+1,nw,1)) / 2.0
     786            endif
    771787            ! Opacity (tau) :
    772788            do k = L_NLAYRAD, lev2lay, -1
    773789               popthv(ig,lev2lay,nw,2) = popthv(ig,lev2lay,nw,2) + popthv(ig,k,nw,1)
    774                poptcv(ig,lev2lay,nw,2) = poptcv(ig,lev2lay,nw,2) + poptcv(ig,k,nw,1)
     790               if (callclouds) then
     791                  poptcv(ig,lev2lay,nw,2) = poptcv(ig,lev2lay,nw,2) + poptcv(ig,k,nw,1)
     792               endif
    775793            enddo
    776794            ! Extinction (k) :
    777795            popthv(ig,lev2lay,nw,3) = popthv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
    778             poptcv(ig,lev2lay,nw,3) = poptcv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
     796            if (callclouds) then
     797               poptcv(ig,lev2lay,nw,3) = poptcv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
     798            endif
    779799            ! Simple Scattering Albedo (w) and Asymmetry Parameter (g) :
    780800            popthv(ig,lev2lay,nw,4) = (popt_hazev(l,nw,2) + popt_hazev(l+1,nw,2)) / 2.0
    781             poptcv(ig,lev2lay,nw,4) = (popt_cloudsv(l,nw,2) + popt_cloudsv(l+1,nw,2)) / 2.0
    782801            popthv(ig,lev2lay,nw,5) = (popt_hazev(l,nw,3) + popt_hazev(l+1,nw,3)) / 2.0
    783             poptcv(ig,lev2lay,nw,5) = (popt_cloudsv(l,nw,3) + popt_cloudsv(l+1,nw,3)) / 2.0         
     802            if (callclouds) then
     803               poptcv(ig,lev2lay,nw,4) = (popt_cloudsv(l,nw,2) + popt_cloudsv(l+1,nw,2)) / 2.0
     804               poptcv(ig,lev2lay,nw,5) = (popt_cloudsv(l,nw,3) + popt_cloudsv(l+1,nw,3)) / 2.0
     805            endif
    784806         enddo
    785 
    786807         ! Infra-Red
    787808         do nw=1,L_NSPECTI
    … …  
    790811            ! Optical thickness (dtau) :
    791812            popthi(ig,lev2lay,nw,1) = (popt_hazei(l,nw,1) + popt_hazei(l+1,nw,1)) / 2.0
    792             poptci(ig,lev2lay,nw,1) = (popt_cloudsi(l,nw,1) + popt_cloudsi(l+1,nw,1)) / 2.0
     813            if (callclouds) then
     814               poptci(ig,lev2lay,nw,1) = (popt_cloudsi(l,nw,1) + popt_cloudsi(l+1,nw,1)) / 2.0
     815            endif
    793816            ! Opacity (tau) :
    794817            do k = L_NLAYRAD, lev2lay, -1
    795818               popthi(ig,lev2lay,nw,2) = popthi(ig,lev2lay,nw,2) + popthi(ig,k,nw,1)
    796                poptci(ig,lev2lay,nw,2) = poptci(ig,lev2lay,nw,2) + poptci(ig,k,nw,1)
     819               if (callclouds) then
     820                  poptci(ig,lev2lay,nw,2) = poptci(ig,lev2lay,nw,2) + poptci(ig,k,nw,1)
     821               endif
    797822            enddo
    798823            ! Extinction (k) :
    799824            popthi(ig,lev2lay,nw,3) = popthi(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
    800             poptci(ig,lev2lay,nw,3) = poptci(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
     825            if (callclouds) then
     826               poptci(ig,lev2lay,nw,3) = poptci(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
     827            endif
    801828            ! Simple Scattering Albedo (w) and Asymmetry Parameter (g) :
    802829            popthi(ig,lev2lay,nw,4) = (popt_hazei(l,nw,2) + popt_hazei(l+1,nw,2)) / 2.0
    803             poptci(ig,lev2lay,nw,4) = (popt_cloudsi(l,nw,2) + popt_cloudsi(l+1,nw,2)) / 2.0
    804830            popthi(ig,lev2lay,nw,5) = (popt_hazei(l,nw,3) + popt_hazei(l+1,nw,3)) / 2.0
    805             poptci(ig,lev2lay,nw,5) = (popt_cloudsi(l,nw,3) + popt_cloudsi(l+1,nw,3)) / 2.0
     831            if (callclouds) then
     832               poptci(ig,lev2lay,nw,4) = (popt_cloudsi(l,nw,2) + popt_cloudsi(l+1,nw,2)) / 2.0
     833               poptci(ig,lev2lay,nw,5) = (popt_cloudsi(l,nw,3) + popt_cloudsi(l+1,nw,3)) / 2.0
     834            endif
    806835         enddo
    807836      enddo
    … …  
    826855                      real(-nfluxtopv),real(nfluxtopi)
    827856         close(116)
    828 
    829857
    830858!          USEFUL COMMENT - Do Not Remove.

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

    r3083 r3090  
    5858      real,save    :: FHIR
    5959!$OMP THREADPRIVATE(Fcloudy,Fssa,FHVIS,FHIR)
     60     
    6061      logical,save :: resCH4_inf
    6162!$OMP THREADPRIVATE(resCH4_inf)

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

    r3083 r3090  
    1818  !! done elsewhere.
    1919  !!
    20   !! Authors : J.Burgalat, J.Vatant d'Ollone - 2017, B. de Batz de Trenquelléon (2023)
     20  !! Authors : J.Burgalat, J.Vatant d'Ollone - 2017
     21  !! Modified : B. de Batz de Trenquelléon (2023)
    2122  !!
    2223  USE MMP_GCM

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

    r3083 r3090  
    2424!   dtlc = Condensation heating rate [K.s-1]
    2525!
    26 !   Author(s)
    27 !   ---------
     26!   Author
     27!   ------
    2828!   B. de Batz de Trenquelléon (07/2023)
    2929!====================================================================   

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

    r3083 r3090  
    1717!     Jan Vatant d'Ollone (2016)
    1818!
     19!  Modified :
     20!     B. de Batz de Trenquelléon (2022)
    1921! ==========================================================================
    2022

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

    r3083 r3090  
    4444!
    4545!
    46 !   Author(s)
    47 !   ---------
     46!   Author
     47!   ------
    4848!   B. de Batz de Trenquelléon (10/2022)
    4949!
    … …  
    7979real, parameter :: humCH4 = 0.5     ! Imposed surface humidity for CH4 [-]
    8080
    81 real, parameter :: Flnp = 0.01      ! Fraction occupied by lakes (North Pole)
     81real, parameter :: Flnp = 0.05      ! Fraction occupied by lakes (North Pole)
    8282real, parameter :: Flsp = 0.005     ! Fraction occupied by lakes (South Pole)
    8383
    … …  
    129129  qsatCH4 = (1.0e5 / pplev(ig,1)) * exp(1.051e1 - 1.110e3/tsurf(ig) - 4.341e3/tsurf(ig)**2 + 1.035e5/tsurf(ig)**3 - 7.910e5/tsurf(ig)**4)
    130130  qsatCH4 = humCH4 * qsatCH4
    131   ! CH4 : 0.85 * qsat because of dissolution in N2
    132   qsatCH4 = 0.85 * qsatCH4
     131  ! CH4 : 0.80 * qsat because of dissolution in N2
     132  qsatCH4 = 0.80 * qsatCH4
    133133
    134134  ! Flux at the surface [kg.m-2.s-1] :
    … …  
    196196  dtsurfevap(ig) = - (fluxCH4 * LvCH4) / Cplake
    197197
    198   ! >>> [TEMPO : BBT]
    199   !open(501,file='Evap_CH4.out')
    200   !if(REAL(latitude_deg(ig)) .ge. 70.) then
    201   !      write(501,*) '-----------------------------------'
    202   !      write(501,*) 'Latitude =', REAL(latitude_deg(ig))
    203   !      write(501,*) ''
    204   !      write(501,*) 'Initialisation variables :'
    205   !     write(501,*) 'Air density =', rhoair, 'kg.m-3'
    206   !      write(501,*) 'Horizontal wind =', ws, 'm.s-1'
    207   !      write(501,*) 'Turbulent term =', Cd
    208   !      write(501,*) 'CH4 saturation at the surface =', qsatCH4, 'mol.mol-1'
    209   !      write(501,*) ''
    210   !      write(501,*) 'Flux :'
    211   !      write(501,*) 'Flux of CH4 =', fluxCH4, 'kg.m-2.s-1.mol.mol-1'
    212   !      write(501,*) ''
    213   !      write(501,*) 'Variation of CH4 :'
    214   !      write(501,*) 'New molar fraction of CH4 =', newpqCH4, 'mol.mol-1'
    215   !      write(501,*) 'Trend of CH4 =', pdqCH4(ig), 'mol.mol-1.s-1'
    216   !      write(501,*) 'Variation of CH4 tank =', dtankCH4, 'm'
    217   !      write(501,*) ''
    218   !      write(501,*) 'Latent Heat :'
    219   !      write(501,*) 'Latente heat of CH4 vaporisation =', LvCH4, 'J.kg-1.mol.mol-1.s-1'
    220   !      write(501,*) 'dtsurfevap :', dtsurfevap(ig), 'K.s-1'
    221   !      write(501,*) '-----------------------------------'
    222   !endif
    223   !close(501)
    224   ! <<< [TEMPO : BBT]
    225 
    226198enddo ! End of main loop on the horizontal grid
    227199

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

    r3083 r3090  
    4343  !
    4444  !
    45   !     Authors
    46   !     -------
     45  !     Author
     46  !     ------
    4747  !     B. de Batz de Trenquelleon (08/2022).
    4848  !

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

    r3083 r3090  
    3333!   author: Frederic Hourdin 15 / 10 /93
    3434!   -------
    35 !   modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
    36 !             Ehouarn Millour (oct. 2008) tracers are now identified
    37 !              by their names and may not be contiguously
    38 !              stored in the q(:,:,:,:) array
    39 !             E.M. (june 2009) use getin routine to load parameters
     35!   modified: -Sebastien Lebonnois 11/06/2003 (new callphys.def)
     36!             -Ehouarn Millour (oct. 2008) tracers are now identified
     37!               by their names and may not be contiguously
     38!               stored in the q(:,:,:,:) array
     39!             -Ehouarn Millour (june 2009) use getin routine
     40!               to load parameters
     41!             -Bruno de Batz de Trenquelléon (2023) minor changes and
     42!               new initialisations
    4043!
    4144!

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

    r3083 r3090  
    2323  !     -------
    2424  !     J. Vatant d'Ollone, J.Burgalat, S.Lebonnois (09/2017)
    25   !     B. de Batz de Trenquelléon (02/2023)
     25  !     Modified : B. de Batz de Trenquelléon (02/2023)
    2626  !
    2727  !============================================================================

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

    r3083 r3090  
    3333  !     -------
    3434  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
    35   !     Clean and adaptation to Titan by J. Vatant d'Ollone (2016-17)
    36   !     Clean and correction to Titan by B. de Batz de Trenquelléon (2022)
    37   !     New optics added for Titan's clouds by B. de Batz de Trenquelléon (2023)
    38   !     The whole routine needs to be redone...
     35  !     
     36  !     Modified
     37  !     --------
     38  !     J. Vatant d'Ollone (2016-17)
     39  !         --> Clean and adaptation to Titan
     40  !     B. de Batz de Trenquelléon (2022-2023)
     41  !         --> Clean and correction to Titan
     42  !         --> New optics added for Titan's clouds
    3943  !     
    4044  !==================================================================
    … …  
    105109  ! Variables for new optics
    106110  integer iq, iw, FTYPE, CTYPE
    107   real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3ices,m3cld
     111  real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3cld
    108112  real*8 dtauaer_s,dtauaer_f,dtau_ccn,dtau_cld
    109113  real*8,save :: rhoaer_s(L_NSPECTI),ssa_s(L_NSPECTI),asf_s(L_NSPECTI)
    … …  
    178182               m3as = pqmo(ilay,2) / 2.0
    179183               m3af = pqmo(ilay,4) / 2.0
    180                ! Cut-off
    181                IF (ilay .lt. 19) THEN
    182                   m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
    183                   m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     184               ! Cut-off (here for p = 2.7e3Pa / alt = 70km)
     185               IF (ilay .lt. 23) THEN
     186                  m3as = pqmo(23,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
     187                  m3af = pqmo(23,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
    184188               ENDIF
    185189
    … …  
    202206               m0as  = pqmo(ilay,1) / 2.0
    203207               m3as  = pqmo(ilay,2) / 2.0
    204                ! If not callclouds : must have a cut-off
     208               ! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km)
    205209               IF (.NOT. callclouds) THEN
    206                   IF (ilay .lt. 19) THEN
    207                      m0as = pqmo(19,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
    208                      m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     210                  IF (ilay .lt. 23) THEN
     211                     m0as = pqmo(23,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
     212                     m3as = pqmo(23,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
    209213                  ENDIF
    210214               ENDIF
    … …  
    216220               m0af  = pqmo(ilay,3) / 2.0
    217221               m3af  = pqmo(ilay,4) / 2.0
    218                ! If not callclouds : must have a cut-off
     222               ! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km)
    219223               IF (.NOT. callclouds) THEN
    220                   IF (ilay .lt. 19) THEN
    221                      m0af = pqmo(19,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
    222                      m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     224                  IF (ilay .lt. 23) THEN
     225                     m0af = pqmo(23,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
     226                     m3af = pqmo(23,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
    223227                  ENDIF
    224228               ENDIF
    … …  
    249253               m0ccn = pqmo(ilay,5) / 2.0
    250254               m3ccn = pqmo(ilay,6) / 2.0
    251                m3ices = 0.0d0
    252                m3cld  = m3ccn
     255               m3cld  = 0.0d0
    253256               
    254257               ! Clear / Dark column method :
    … …  
    261264               IF (CDCOLUMN == 0) THEN
    262265                  DO iq = 2, nice
    263                      m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
    264266                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
    265267                  ENDDO
    … …  
    269271               ELSEIF (CDCOLUMN == 1) THEN
    270272                  DO iq = 1, nice
    271                      m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
    272273                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
    273274                  ENDDO
    … …  
    278279                  WRITE(*,*) 'WE USE DARK COLUMN ...'
    279280                  DO iq = 1, nice
    280                      m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
    281281                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
    282282                  ENDDO
    … …  
    285285               
    286286               ! For small dropplets, opacity of nucleus dominates...
    287                !ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3ices) / (m3ccn + m3ices)
    288287               ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3cld) / (m3ccn + m3cld)
    289288

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

    r3083 r3090  
    2525  !     -------
    2626  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
    27   !     Clean and adaptation to Titan by J. Vatant d'Ollone (2016-17)
    28   !     Clean and correction to Titan by B. de Batz de Trenquelléon (2022)
    29   !     New optics added for Titan's clouds by B. de Batz de Trenquelléon (2023)
    30   !     The whole routine needs to be redone...
     27  !     
     28  !     Modified
     29  !     --------
     30  !     J. Vatant d'Ollone (2016-17)
     31  !         --> Clean and adaptation to Titan
     32  !     B. de Batz de Trenquelléon (2022-2023)
     33  !         --> Clean and correction to Titan
     34  !         --> New optics added for Titan's clouds
    3135  !     
    3236  !==================================================================
    … …  
    116120  ! Variables for new optics
    117121  integer iq, iw, FTYPE, CTYPE
    118   real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3ices,m3cld
     122  real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3cld
    119123  real*8 dtauaer_s,dtauaer_f,dtau_ccn,dtau_cld
    120124  real*8,save :: rhoaer_s(L_NSPECTV),ssa_s(L_NSPECTV),asf_s(L_NSPECTV)
    … …  
    202206               m3as = pqmo(ilay,2) / 2.0
    203207               m3af = pqmo(ilay,4) / 2.0
    204                ! Cut-off
    205                IF (ilay .lt. 19) THEN
    206                   m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
    207                   m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     208               ! Cut-off (here for p = 2.7e3Pa / alt = 70km)
     209               IF (ilay .lt. 23) THEN
     210                  m3as = pqmo(23,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
     211                  m3af = pqmo(23,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
    208212               ENDIF
    209213               
    … …  
    226230               m0as  = pqmo(ilay,1) / 2.0
    227231               m3as  = pqmo(ilay,2) / 2.0
    228                ! If not callclouds : must have a cut-off
     232               ! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km)
    229233               IF (.NOT. callclouds) THEN
    230                   IF (ilay .lt. 19) THEN
    231                      m0as = pqmo(19,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
    232                      m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     234                  IF (ilay .lt. 23) THEN
     235                     m0as = pqmo(23,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
     236                     m3as = pqmo(23,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
    233237                  ENDIF
    234238               ENDIF
    … …  
    240244               m0af  = pqmo(ilay,3) / 2.0
    241245               m3af  = pqmo(ilay,4) / 2.0
    242                ! If not callclouds : must have a cut-off
     246               ! If not callclouds : must have a cut-off (here for p = 2.7e3Pa / alt = 70km)
    243247               IF (.NOT. callclouds) THEN
    244                   IF (ilay .lt. 19) THEN
    245                      m0af = pqmo(19,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
    246                      m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     248                  IF (ilay .lt. 23) THEN
     249                     m0af = pqmo(23,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
     250                     m3af = pqmo(23,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(24)-zlev(23))
    247251                  ENDIF
    248252               ENDIF
    … …  
    273277               m0ccn = pqmo(ilay,5) / 2.0
    274278               m3ccn = pqmo(ilay,6) / 2.0
    275                m3ices = 0.0d0
    276                m3cld  = m3ccn
     279               m3cld = 0.0d0
    277280               
    278281               ! Clear / Dark column method :
    … …  
    285288               IF (CDCOLUMN == 0) THEN
    286289                  DO iq = 2, nice
    287                      m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
    288290                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
    289291                  ENDDO
    … …  
    293295               ELSEIF (CDCOLUMN == 1) THEN
    294296                  DO iq = 1, nice
    295                      m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
    296297                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
    297298                  ENDDO
    … …  
    302303                  WRITE(*,*) 'WE USE DARK COLUMN ...'
    303304                  DO iq = 1, nice
    304                      m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
    305305                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
    306306                  ENDDO
    … …  
    309309
    310310               ! For small dropplets, opacity of nucleus dominates
    311                !ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3ices) / (m3ccn + m3ices)
    312311               ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3cld) / (m3ccn + m3cld)
    313312               ssa_cld(nw) = ssa_cld(nw) * Fssa

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

    r3083 r3090  
    1515 
    1616      use radinc_h, only : L_NSPECTI,L_NSPECTV
    17       use radcommon_h, only: sigma, gzlat, grav, BWNV
     17      use radcommon_h, only: sigma, gzlat, grav, BWNV, WNOI, WNOV
    1818      use comchem_h, only: nkim, cnames, nlaykim_up, ykim_up, ykim_tot, botCH4
    1919      use comdiurn_h, only: coslat, sinlat, coslon, sinlon
    … …  
    161161!                + clean of all too-generic (ocean, water, co2 ...) routines
    162162!                + Titan's chemistry
    163 !           Microphysical moment model - J.Burgalat / J.Vatant d'Ollone (2017-2018)
    164 !           Improvement of microphysical moment model - B. de Batz de Trenquelléon (2022-2023)
    165 !           Optics ofnn Titan's coulds - B. de Batz de Trenquelléon (2022-2023)
     163!           Microphysical moment model - J.Burgalat / B. de Batz de Trenquelléon (2022-2023)
     164!           Modified : B. de Batz de Trenquelléon (2023)
     165!                + Optical improvements (haze and cloud)
     166!                + Nudging of troposphere !!
    166167!============================================================================================
    167168
    … …  
    258259      real zdqmr(ngrid,nlayer,nq)     ! Mass_redistribution routine.
    259260     
    260       real zdqchi(ngrid,nlayer,nq)    ! Chemical tendency ( chemistry routine ).
     261      real zdqchi(ngrid,nlayer,nq)    ! Chemical tendency (chemistry routine).
    261262     
    262263      real zdqmufi(ngrid,nlayer,nq)   ! Microphysical tendency.
    … …  
    271272      real zdhdif(ngrid,nlayer)                      ! Turbdiff/vdifc routines.
    272273      real zdhadj(ngrid,nlayer)                      ! Convadj routine.
    273       real zdundg(ngrid,nlayer)                      ! Nudging for zonal wind (m/s/s).
     274      real zdundg(ngrid,nlayer)                      ! Nudging for zonal wind.
    274275     
    275276      ! For Pressure and Mass :
    … …  
    301302      real zdudyn(ngrid,nlayer)                          ! Dynamical Zonal Wind tendency (m.s-2).
    302303
    303       real zhorizwind(ngrid,nlayer) ! Horizontal Wind ( sqrt(u*u+v*v))
     304      real zhorizwind(ngrid,nlayer) ! Horizontal Wind (sqrt(u*u+v*v))
    304305
    305306      real vmr(ngrid,nlayer)                        ! volume mixing ratio
    … …  
    352353      real, dimension(ngrid,nlayer,nkim) :: ysat
    353354     
    354       ! Fraction of CH4 [mol/mol]
     355      ! Temporary fraction of CH4 [mol/mol]
    355356      real, dimension(ngrid,nlayer) :: tpq_CH4
    356357
    … …  
    467468!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    468469
    469          if ( callchim ) then
    470 
    471             if ( moyzon_ch .and. ngrid.eq.1 ) then
     470         if (callchim) then
     471
     472            if (moyzon_ch .and. ngrid.eq.1) then
    472473               print *, "moyzon_ch=",moyzon_ch," and ngrid=1"
    473474               print *, "Please desactivate zonal mean for 1D !"
    … …  
    483484!        ~~~~~~~~~~~~~~~~~~~~~~~~
    484485
    485          IF ( callmufi ) THEN
     486         IF (callmufi) THEN
    486487           ! WARNING JB (27/03/2018): inimufi should be wrapped in an OMP SINGLE statement.
    487488           call inimufi(ptimestep)
    … …  
    619620         write(*,*) "physiq: call initialize_xios_output"
    620621         call initialize_xios_output(pday,ptime,ptimestep,daysec, &
    621                                      presnivs,pseudoalt)
     622                                     presnivs,pseudoalt,wnoi,wnov)
    622623#endif
    623624         write(*,*) "physiq: end of firstcall"
    … …  
    10661067      if (tracer) then
    10671068
    1068 
    10691069#ifndef MESOSCALE
    1070 !! JVO 20 : For now, no chemistry or microphysics in MESOSCALE, but why not in the future ?
    1071 
    1072   ! -------------------
    1073   !   V.1 Microphysics
    1074   ! -------------------
    1075 
    1076          ! JVO 05/18 : We must call microphysics before chemistry, for condensation !
    1077  
     1070   ! -------------------
     1071   !   V.1 Microphysics
     1072   ! -------------------
     1073
     1074         ! We must call microphysics before chemistry, for condensation !
    10781075         if (callmufi) then
    10791076
    … …  
    10841081            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,tpq,dtpq,zdqmufi)
    10851082            tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)*ptimestep ! only manipulation of tpq->*ptimestep here
    1086 #else
    1087             
     1083
     1084#else   
    10881085            call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi)
    1089 
    10901086            pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
    10911087
    1092             ! Sanity check ( way safer to be done here rather than within YAMMS )
     1088            ! Sanity check (way safer to be done here rather than within YAMMS)
    10931089            ! Important : the sanity check intentionally include the former processes tendency !
    10941090            ! NB : Despite this sanity check there might be still some unphysical values going through :
    … …  
    11241120
    11251121            ! Microphysics condensation for 2D fields to sent non-saturated fields to photochem
    1126             if ( callclouds .and. moyzon_ch .and. mod(icount-1,ichim).eq.0 ) then
    1127               zdqfibar(:,:,:) = 0.D0 ! We work in zonal average -> forget processes other than condensation
    1128               call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, &
    1129                            gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar)
    1130               ! TODO : Add a sanity check here !
     1122            if (callclouds .and. moyzon_ch .and. mod(icount-1,ichim).eq.0) then
     1123               zdqfibar(:,:,:) = 0.D0 ! We work in zonal average -> forget processes other than condensation
     1124               call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, &
     1125                            gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar)
     1126               ! TODO : Add a sanity check here !
    11311127            endif
    11321128         
    1133             ! >>> TEMPO : BBT
    1134             ! Default value -> no condensation [kg/kg_air/s] :
    1135             dmuficond(:,:,:) = 0.D0
    1136             do iq = 1, size(ices_indx)
    1137                dmuficond(:,:,iq) = zdqmufi(:,:,gazs_indx(iq))
    1138             enddo
    1139             call cond_muphy(ngrid,nlayer,pt,dmuficond,zdtlc)
    1140             !pdt(:,:) = pdt(:,:) + zdtlc(:,:)
    1141             ! <<< TEMPO : BBT
    1142          endif
     1129            ! Condensation heating rate :
     1130            if (callclouds) then
     1131               ! Default value -> no condensation [kg/kg_air/s] :
     1132               dmuficond(:,:,:) = 0.D0
     1133               do iq = 1, size(ices_indx)
     1134                  dmuficond(:,:,iq) = zdqmufi(:,:,gazs_indx(iq))
     1135               enddo
     1136               call cond_muphy(ngrid,nlayer,pt,dmuficond,zdtlc)
     1137               !pdt(:,:) = pdt(:,:) + zdtlc(:,:)
     1138            endif
     1139         endif ! callmufi
    11431140     
    11441141  ! -----------------
    … …  
    11461143  ! -----------------
    11471144  !   NB : Must be call last ( brings fields back to an equilibrium )
    1148 
    1149   ! Known bug ? ( JVO 18 ) : If you try to use chimi_indx instead of iq+nmicro
    1150   ! it leads to weird results / crash on dev mod ( ok in debug ) ... Why ? Idk ...
    1151 
    11521145         if (callchim) then
    11531146
    … …  
    11551148            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    11561149            do iq = 1,nkim
    1157                ychim(:,:,iq) = ( pq(:,:,iq+nmicro) + pdq(:,:,iq+nmicro)*ptimestep ) / rat_mmol(iq+nmicro)
     1150               ychim(:,:,iq) = (pq(:,:,iq+nmicro) + pdq(:,:,iq+nmicro)*ptimestep) / rat_mmol(iq+nmicro)
    11581151            enddo
    11591152
    … …  
    11621155            ! processes ( convadj ... ) we miss them in any case as we work in zonally/diurnal
    11631156            ! mean -> no fine diurnal/short time evolution, only seasonal evolution only.
    1164             if ( moyzon_ch .and. mod(icount-1,ichim).eq. 0 ) then
     1157            if (moyzon_ch .and. mod(icount-1,ichim).eq. 0) then
    11651158              do iq = 1,nkim
    11661159                ychimbar(:,:,iq) =  zqfibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
    … …  
    11741167            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    11751168           
    1176             call calc_ysat(ngrid,nlayer,pplay/100.0,pt,ysat) ! Compute saturation profiles for every grid point
     1169            call calc_ysat(ngrid,nlayer,pplay/100.0,pt,ysat) ! Compute saturation profiles for every grid point (!!p in mbar!!)
    11771170
    11781171            dyccond(:,:,:) = 0.D0 ! Default value -> no condensation
    11791172
    11801173            do iq=1,nkim
    1181                where ( ychim(:,:,iq).gt.ysat(:,:,iq) )   &
    1182                      dyccond(:,:,iq+nmicro) = ( -ychim(:,:,iq)+ysat(:,:,iq) ) / ptimestep
     1174               where (ychim(:,:,iq).gt.ysat(:,:,iq))
     1175                  dyccond(:,:,iq+nmicro) = (-ychim(:,:,iq)+ysat(:,:,iq)) / ptimestep
     1176               endwhere
    11831177            enddo
    11841178
    … …  
    11911185            do iq=1,nkim
    11921186              ychim(:,:,iq) = ychim(:,:,iq) + dyccond(:,:,iq+nmicro)*ptimestep ! update molar ychim for following calchim
    1193 
    11941187              pdq(:,:,iq+nmicro) = pdq(:,:,iq+nmicro) + dyccond(:,:,iq+nmicro)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
    11951188            enddo
    … …  
    11981191            ! ii. 2D zonally averaged fields need to condense and evap before photochem
    11991192            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1200             if ( moyzon_ch .and. mod(icount-1,ichim).eq. 0 ) then
    1201 
    1202               call calc_ysat(ngrid,nlayer,zplaybar/100.0,ztfibar,ysat) ! Compute saturation profiles for every grid point for the zon-ave fields
    1203 
    1204               dyccondbar(:,:,:) = 0.D0 ! Default value -> no condensation
    1205              
    1206               do iq = 1,nkim
    1207                  where ( ychimbar(:,:,iq).gt.ysat(:,:,iq) )  &
    1208                        dyccondbar(:,:,iq+nmicro) = ( -ychimbar(:,:,iq)+ysat(:,:,iq) ) / ptimestep
    1209               enddo
     1193            if (moyzon_ch .and. mod(icount-1,ichim).eq. 0) then
     1194
     1195               call calc_ysat(ngrid,nlayer,zplaybar/100.0,ztfibar,ysat) ! Compute saturation profiles for every grid point for the zon-ave fields
     1196
     1197               dyccondbar(:,:,:) = 0.D0 ! Default value -> no condensation
     1198               
     1199               do iq = 1,nkim
     1200                     where ( ychimbar(:,:,iq).gt.ysat(:,:,iq) )
     1201                        dyccondbar(:,:,iq+nmicro) = ( -ychimbar(:,:,iq)+ysat(:,:,iq) ) / ptimestep
     1202                     endwhere
     1203               enddo
    12101204
    12111205               if (callclouds) then
    … …  
    12151209               endif
    12161210
    1217               do iq=1,nkim
    1218                 ychimbar(:,:,iq) = ychimbar(:,:,iq) + dyccondbar(:,:,iq+nmicro)*ptimestep
    1219               enddo
    1220 
    1221               ! Pseudo-evap ( forcing constant surface humidity )
    1222               !do ig=1,ngrid
    1223               !   if ( ychimbar(ig,1,7) .lt. botCH4 ) ychimbar(ig,1,7) = botCH4
    1224               !enddo
    1225 
    1226             endif ! if ( moyzon_ch .and. mod(icount-1,ichim).eq. 0 )
     1211               do iq=1,nkim
     1212                  ychimbar(:,:,iq) = ychimbar(:,:,iq) + dyccondbar(:,:,iq+nmicro)*ptimestep
     1213               enddo
     1214
     1215            endif ! if ( moyzon_ch .and. mod(icount-1,ichim).eq.0 )
    12271216
    12281217            ! iii. Photochemistry ( must be call after condensation (and evap of 2D) )
    12291218            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    12301219            if( mod(icount-1,ichim).eq.0. ) then
    1231                
    1232                print *, "We enter in the chemistry ..."
     1220               print *, "We enter in the photochemistry ..."
    12331221
    12341222               if (moyzon_ch) then ! 2D zonally averaged chemistry
    1235 
    12361223                 ! Here we send zonal average fields ( corrected with cond ) from dynamics to chem. module
    12371224                 call calchim(ngrid,ychimbar,declin,ctimestep,ztfibar,zphibar,zphisbar,  &
    12381225                              zplaybar,zplevbar,zzlaybar,zzlevbar,dycchi)
     1226               
    12391227               else ! 3D chemistry (or 1D run)
    12401228                 call calchim(ngrid,ychim,declin,ctimestep,pt,pphi,phisfi,  &
    12411229                              pplay,pplev,zzlay_eff,zzlev_eff,dycchi)
    12421230               endif ! if moyzon
    1243 
    1244             endif
     1231            endif ! if (mod(icount-1,ichim).eq.0)
    12451232           
    12461233            ! iv. Surface pseudo-evaporation
    12471234            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     1235            ! Infinite tank of CH4 (global)
    12481236            if (resCH4_inf) then
    1249                 do ig=1,ngrid
    1250                         if ( (ychim(ig,1,7)+dycchi(ig,1,7)*ptimestep) .lt. botCH4 ) then ! + dycchi because ychim not yet updated
    1251                                 dycevapCH4(ig) = ( -ychim(ig,1,7)+botCH4 ) / ptimestep - dycchi(ig,1,7)
    1252                         else
    1253                                 dycevapCH4(ig) = 0.D0
    1254                         endif
    1255                 enddo
    1256 
    1257             ! >>> New evaporation flux
     1237               do ig=1,ngrid
     1238                  if ((ychim(ig,1,7)+dycchi(ig,1,7)*ptimestep) .lt. botCH4) then ! + dycchi because ychim not yet updated
     1239                     dycevapCH4(ig) = ( -ychim(ig,1,7)+botCH4 ) / ptimestep - dycchi(ig,1,7)
     1240                  else
     1241                     dycevapCH4(ig) = 0.D0
     1242                  endif
     1243               enddo
     1244
     1245            ! Tank at the pole (for the moment infinit...)
    12581246            else
    1259                 tankCH4(:) = 2. ! Pour le moment reservoir infini...
    1260                 if (moyzon_ch) then
    1261                     tpq_CH4(:,:) = ychimbar(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
    1262                 else
    1263                     tpq_CH4(:,:) = ychim(:,:,7) + dycchi(:,:,7)*ptimestep    ! + dycchi because ychim not yet updated [mol/mol]
    1264                 endif
    1265                 call evapCH4(ngrid,nlayer,ptimestep,pplev,zzlay,zzlev,&
    1266                              pu,pv,tsurf,tpq_CH4,tankCH4,dycevapCH4,zdtsurfevap)
    1267                 zdtsurf(:) = zdtsurf(:) + zdtsurfevap(:)
     1247               tankCH4(:) = 2.
     1248               if (moyzon_ch) then
     1249                 tpq_CH4(:,:) = ychimbar(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
     1250               else
     1251                 tpq_CH4(:,:) = ychim(:,:,7) + dycchi(:,:,7)*ptimestep    ! + dycchi because ychim not yet updated [mol/mol]
     1252               endif
     1253               call evapCH4(ngrid,nlayer,ptimestep,pplev,zzlay,zzlev,&
     1254                            pu,pv,tsurf,tpq_CH4,tankCH4,dycevapCH4,zdtsurfevap)
     1255               zdtsurf(:) = zdtsurf(:) + zdtsurfevap(:)
    12681256            endif
    1269             ! <<< New evaporation flux
    12701257
    12711258            pdq(:,1,7+nmicro) = pdq(:,1,7+nmicro) + dycevapCH4(:)*rat_mmol(7+nmicro)
    … …  
    12861273         endif ! end of 'callchim'
    12871274
     1275! END MESOSCALE
    12881276#endif
    12891277
    … …  
    13161304            pdv(:,:)    = pdv(:,:)    + zdvmr(:,:)
    13171305            pdpsrf(:)   = pdpsrf(:)   + zdpsrfmr(:)
     1306            zdtsurf(:)  = zdtsurf(:)  + zdtsurfmr(:)
    13181307         endif
    13191308
    … …  
    13621351      if (nudging_u) then
    13631352         zdundg(:,:) = 0.D0
    1364          
    13651353         ! boucle sur les points de grille :
    13661354         do i = 1, ngrid
    13671355            ! boucle sur les latitudes :
    13681356            do j = 1, 49
    1369 
    13701357               ! Nudging of the first 23 layers only !!
    13711358               if (ABS(REAL(u_ref(j,1)) - REAL(latitude_deg(i))) .lt. 1e-2) then
    13721359                  zdundg(i,1:23) = (u_ref(j,2:24) - (pu(i,1:23)+pdu(i,1:23)*ptimestep)) / nudging_dt
    13731360               endif
    1374            
    13751361            enddo
    13761362         enddo
    1377 
    13781363         pdu(:,:) = pdu(:,:) + zdundg(:,:)
    13791364      endif
    … …  
    14061391      enddo
    14071392
    1408       ! Tracers.
    1409       ! >>> TEMPO : BBT
    1410       pdq(:,1:13,gazs_indx(2)) = (2e-8 - pq(:,1:13,gazs_indx(2))) / (ptimestep * 100)  ! C2H6
    1411       pdq(:,1:12,gazs_indx(3)) = (1e-9 - pq(:,1:12,gazs_indx(3))) / (ptimestep * 100)  ! C2H2
    1412       pdq(:,1:11,gazs_indx(4)) = (1e-20 - pq(:,1:11,gazs_indx(4))) / (ptimestep * 100) ! HCN
    1413       ! <<< TEMPO : BBT
     1393      ! Tracers : for now because problem with re-evaporation in microphysics ...
     1394      if (callclouds) then
     1395         pdq(:,1:13,gazs_indx(2)) = (2e-8 - pq(:,1:13,gazs_indx(2))) / (ptimestep * 100)  ! C2H6
     1396         pdq(:,1:12,gazs_indx(3)) = (1e-9 - pq(:,1:12,gazs_indx(3))) / (ptimestep * 100)  ! C2H2
     1397         pdq(:,1:11,gazs_indx(4)) = (1e-20 - pq(:,1:11,gazs_indx(4))) / (ptimestep * 100) ! HCN
     1398      endif
    14141399      zq(:,:,:) = pq(:,:,:) + pdq(:,:,:)*ptimestep
    14151400
    … …  
    15011486      endif
    15021487
    1503 
    15041488      if (is_master) print*,'--> Ls =',zls*180./pi
    15051489     
    … …  
    15551539      call writediagfi(ngrid,"w","Vertical wind","m.s-1",3,pw)
    15561540      call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
    1557 
    1558       ! Total energy balance diagnostics
     1541     
     1542      ! Subsurface temperatures
     1543      !call writediagsoil(ngrid,"tempsoil","temperature soil","K",3,tsoil)
     1544
     1545     ! Total energy balance diagnostics
    15591546      if(callrad.and.(.not.newtonian))then
    15601547         call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
    … …  
    16661653      ! Send fields to XIOS: (NB these fields must also be defined as
    16671654      ! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
    1668      
     1655
    16691656      !--------------------------------------------------------
    16701657      ! General diagnostics :
    … …  
    16971684      ENDIF
    16981685
     1686      ! Subsurface temperatures (3D)
     1687      !CALL send_xios_field("tempsoil",tsoil)
     1688      ! Total energy balance diagnostics (2D)
     1689      !CALL send_xios_field("ALB",albedo_equivalent)
    16991690
    17001691      !--------------------------------------------------------
    … …  
    17091700      CALL send_xios_field("dudif",zdudif)
    17101701      CALL send_xios_field("duadj",zduadj)
    1711       CALL send_xios_field("dundg",zdundg)
     1702      IF (nudging_u) THEN
     1703         CALL send_xios_field("dundg",zdundg)
     1704      ENDIF
    17121705     
    17131706      ! zhorizwind = sqrt(u*u + v*v)
    … …  
    17271720      CALL send_xios_field("dtdif",zdtdif)
    17281721      CALL send_xios_field("dtadj",zdtadj(:,:))
    1729       CALL send_xios_field("dtlc",zdtlc)
     1722      IF (callclouds) THEN
     1723         CALL send_xios_field("dtlc",zdtlc)
     1724      ENDIF
    17301725
    17311726      ! dtrad = zdtsw + zdtlw
    … …  
    17521747      ! Optical diagnostics :
    17531748      !--------------------------------------------------------
    1754       ! Optics diagnostics for haze (Channel VI04 : 2.80um).
    1755       CALL send_xios_field('tauhv_04',zpopthv(:,:,4,2))
    1756       CALL send_xios_field('khv_04',zpopthv(:,:,4,3))
    1757       CALL send_xios_field('wbarhv_04',zpopthv(:,:,4,4))
    1758       CALL send_xios_field('gbarhv_04',zpopthv(:,:,4,5))
    1759       ! Optics diagnostics for haze (Channel VI18 : 0.85um).
    1760       CALL send_xios_field('tauhv_18',zpopthv(:,:,18,2))
    1761       CALL send_xios_field('khv_18',zpopthv(:,:,18,3))
    1762       CALL send_xios_field('wbarhv_18',zpopthv(:,:,18,4))
    1763       CALL send_xios_field('gbarhv_18',zpopthv(:,:,18,5))
    1764       ! Optics diagnostics for haze (Channel IR07 : 29.5um).
    1765       CALL send_xios_field('tauhi_07',zpopthi(:,:,7,2))
    1766       CALL send_xios_field('khi_07',zpopthi(:,:,7,3))
    1767       CALL send_xios_field('wbarhi_07',zpopthi(:,:,7,4))
    1768       CALL send_xios_field('gbarhi_07',zpopthi(:,:,7,5))
    1769       ! Optics diagnostics for haze (Channel IR22 : 5.66um).
    1770       CALL send_xios_field('tauhi_22',zpopthi(:,:,22,2))
    1771       CALL send_xios_field('khi_22',zpopthi(:,:,22,3))
    1772       CALL send_xios_field('wbarhi_22',zpopthi(:,:,22,4))
    1773       CALL send_xios_field('gbarhi_22',zpopthi(:,:,22,5))
    1774 
    1775       ! Optics diagnostics for haze + clouds (Channel VI04 : 2.80um).
    1776       CALL send_xios_field('taucv_04',zpoptcv(:,:,4,2))
    1777       CALL send_xios_field('kcv_04',zpoptcv(:,:,4,3))
    1778       CALL send_xios_field('wbarcv_04',zpoptcv(:,:,4,4))
    1779       CALL send_xios_field('gbarcv_04',zpoptcv(:,:,4,5))
    1780       ! Optics diagnostics for haze + clouds (Channel VI18 : 0.85um).
    1781       CALL send_xios_field('taucv_18',zpoptcv(:,:,18,2))
    1782       CALL send_xios_field('kcv_18',zpoptcv(:,:,18,3))
    1783       CALL send_xios_field('wbarcv_18',zpoptcv(:,:,18,4))
    1784       CALL send_xios_field('gbarcv_18',zpoptcv(:,:,18,5))
    1785       ! Optics diagnostics for haze + clouds (Channel IR07 : 29.5um).
    1786       CALL send_xios_field('tauci_07',zpoptci(:,:,7,2))
    1787       CALL send_xios_field('kci_07',zpoptci(:,:,7,3))
    1788       CALL send_xios_field('wbarci_07',zpoptci(:,:,7,4))
    1789       CALL send_xios_field('gbarci_07',zpoptci(:,:,7,5))
    1790       ! Optics diagnostics for haze + clouds (Channel IR22 : 5.66um).
    1791       CALL send_xios_field('tauci_22',zpoptci(:,:,22,2))
    1792       CALL send_xios_field('kci_22',zpoptci(:,:,22,3))
    1793       CALL send_xios_field('wbarci_22',zpoptci(:,:,22,4))
    1794       CALL send_xios_field('gbarci_22',zpoptci(:,:,22,5))
    1795 
    1796       ! >>> [TEMPO : BBT]
    1797       !DO nw = 1, L_NSPECTV
    1798       !   WRITE(str2,'(i2.2)') nw
    1799       !   CALL send_xios_field('khv'//TRIM(nw),zpopthv(:,:,nw,3))
    1800       !ENDDO
    1801       ! <<< [TEMPO : BBT]
    1802      
     1749      ! Diagnostics for haze :
     1750      ! IR
     1751      CALL send_xios_field('dtauhi',zpopthi(:,:,:,1))
     1752      CALL send_xios_field('tauhi',zpopthi(:,:,:,2))
     1753      CALL send_xios_field('khi',zpopthi(:,:,:,3))
     1754      CALL send_xios_field('whi',zpopthi(:,:,:,4))
     1755      CALL send_xios_field('ghi',zpopthi(:,:,:,5))
     1756      ! VI
     1757      CALL send_xios_field('dtauhv',zpopthv(:,:,:,1))
     1758      CALL send_xios_field('tauhv',zpopthv(:,:,:,2))
     1759      CALL send_xios_field('khv',zpopthv(:,:,:,3))
     1760      CALL send_xios_field('whv',zpopthv(:,:,:,4))
     1761      CALL send_xios_field('ghv',zpopthv(:,:,:,5))
     1762
     1763      ! Diagnostics for haze and clouds :
     1764      IF (callclouds) THEN
     1765         ! IR
     1766         CALL send_xios_field('dtaui',zpoptci(:,:,:,1))
     1767         CALL send_xios_field('taui',zpoptci(:,:,:,2))
     1768         CALL send_xios_field('ki',zpoptci(:,:,:,3))
     1769         CALL send_xios_field('wi',zpoptci(:,:,:,4))
     1770         CALL send_xios_field('gi',zpoptci(:,:,:,5))
     1771         ! VI
     1772         CALL send_xios_field('dtauv',zpoptcv(:,:,:,1))
     1773         CALL send_xios_field('tauv',zpoptcv(:,:,:,2))
     1774         CALL send_xios_field('kv',zpoptcv(:,:,:,3))
     1775         CALL send_xios_field('wv',zpoptcv(:,:,:,4))
     1776         CALL send_xios_field('gv',zpoptcv(:,:,:,5))
     1777      ENDIF     
    18031778
    18041779      !--------------------------------------------------------
    … …  
    18311806            CALL send_xios_field("vsed_ccn",mmd_ccn_w(:,:))
    18321807            CALL send_xios_field("flux_ccn",mmd_ccn_flux(:,:))
    1833             CALL send_xios_field("flux_iCH4",mmd_ice_fluxes(:,:,1))
    1834             CALL send_xios_field("flux_iC2H2",mmd_ice_fluxes(:,:,2))
    1835             CALL send_xios_field("flux_iC2H6",mmd_ice_fluxes(:,:,3))
    1836             CALL send_xios_field("flux_iHCN",mmd_ice_fluxes(:,:,3))
    18371808            DO iq = 1, size(ices_indx)
     1809               CALL send_xios_field('flux_i'//TRIM(nameOfTracer(gazs_indx(iq))),mmd_ice_fluxes(:,:,iq))
    18381810               CALL send_xios_field(TRIM(nameOfTracer(gazs_indx(iq)))//'_sat',mmd_gazs_sat(:,:,iq))
    18391811            ENDDO
    … …  
    18431815         CALL send_xios_field("aer_prec",mmd_aer_prec(:))
    18441816         CALL send_xios_field("ccn_prec",mmd_ccn_prec(:))
    1845          CALL send_xios_field("iCH4_prec",mmd_ice_prec(:,1))
    1846          CALL send_xios_field("iC2H2_prec",mmd_ice_prec(:,2))
    1847          CALL send_xios_field("iC2H6_prec",mmd_ice_prec(:,3))
    1848          CALL send_xios_field("iHCN_prec",mmd_ice_prec(:,3))
     1817         IF (callclouds) THEN
     1818            DO iq = 1, size(ices_indx)
     1819               CALL send_xios_field('i'//TRIM(nameOfTracer(gazs_indx(iq)))//'_prec',mmd_ice_prec(:,iq))
     1820            ENDDO
     1821         ENDIF
    18491822      ENDIF ! of 'if callmufi'
    18501823
    … …  
    18561829         ! Chemical species :
    18571830         DO iq = 1, nkim
     1831            ! If no cloud : gzs_indx uninitialized
    18581832            CALL send_xios_field(trim(cnames(iq)),zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro)) ! kg/kg -> mol/mol
    18591833         ENDDO
    18601834
    18611835         ! Condensation tendencies from microphysics (mol/mol/s) :
    1862          DO iq = 1, size(ices_indx)
    1863             CALL send_xios_field('dmuficond_'//trim(nameOfTracer(gazs_indx(iq))),dmuficond(:,:,iq)/rat_mmol(gazs_indx(iq))) ! kg/kg/s -> mol/mol/s
    1864          ENDDO
     1836         IF (callclouds) THEN
     1837            DO iq = 1, size(ices_indx)
     1838               CALL send_xios_field('dmuficond_'//trim(nameOfTracer(gazs_indx(iq))),dmuficond(:,:,iq)/rat_mmol(gazs_indx(iq))) ! kg/kg/s -> mol/mol/s
     1839            ENDDO
     1840         ENDIF
    18651841         
    18661842         ! Condensation tendencies (mol/mol/s) :

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

    r3083 r3090  
    4848!     -------
    4949!     Adapted from the NASA Ames code by R. Wordsworth (2009)
    50 !     Last update : B. de Batz de Trenquelleon (2023)
     50!     Last update : B. de Batz de Trenquelléon (2023)
    5151!     
    5252!==================================================================

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

    r3083 r3090  
    1 ! AUTHOR : J. Burgalat (27/03/2018)
    2 ! NOTES : Add gazs_indx by B. de Batz de Trenquelléon (31/05/2022)
     1! AUTHORS : J. Burgalat (2018), B. de Batz de Trenquelléon (2022)
    32!
    43! OpenMP directives in this module are inconsistent:

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

    r1943 r3090  
    1111
    1212 INTERFACE send_xios_field
    13     MODULE PROCEDURE histwrite0d_xios,histwrite2d_xios,histwrite3d_xios
     13    MODULE PROCEDURE histwrite0d_xios,histwrite2d_xios,histwrite3d_xios,histwrite4d_xios
    1414 END INTERFACE
    1515 
    … …  
    1818
    1919  SUBROUTINE initialize_xios_output(day,timeofday,dtphys,daysec,&
    20                                     presnivs,pseudoalt)
     20                                    presnivs,pseudoalt,wnoi,wnov)
    2121!  USE mod_phys_lmdz_para, only: gather, bcast, &
    2222!                                jj_nb, jj_begin, jj_end, ii_begin, ii_end, &
    … …  
    4646  REAL,INTENT(IN) :: presnivs(:) ! vertical grid approximate pressure (Pa)
    4747  REAL,INTENT(IN) :: pseudoalt(:) ! vertical grid approximate altitude (km)
    48  
     48  REAL,INTENT(IN) :: wnoi(:) ! Array of wavelength channels for the infrared.
     49  REAL,INTENT(IN) :: wnov(:) ! Array of wavelength channels for the visible.
    4950 
    5051  INTEGER :: i
    … …  
    7273    CALL xios_set_axis_attr("preskim_tot", n_glo=size(preskim_tot), value=preskim_tot,&
    7374                            unit="Pa",positive="down")
    74                            
     75   
     76    IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for IR_Wavelength"
     77    CALL xios_set_axis_attr("IR_Wavelength", n_glo=size(wnoi), value=wnoi,unit="m",positive="up")
     78
     79    IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for VI_Wavelength"
     80    CALL xios_set_axis_attr("VI_Wavelength", n_glo=size(wnov), value=wnov,unit="m",positive="up")
     81                 
    7582    ! Calculation of pseudo-altitudes is still to be done
    7683    !IF (prt_level>=10) WRITE(lunout,*) "initialize_xios_output: call xios_set_axis_attr for zlaykim"
    … …  
    291298  END SUBROUTINE histwrite3d_xios
    292299
     300!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     301
     302  SUBROUTINE histwrite4d_xios(field_name, field)
     303    USE dimphy, only: klon, klev
     304    USE mod_phys_lmdz_para, only: gather_omp, grid1Dto2D_mpi, &
     305                                  jj_nb, klon_mpi
     306    USE xios, only: xios_send_field
     307    USE print_control_mod, ONLY: prt_level,lunout
     308    USE mod_grid_phy_lmdz, ONLY: nbp_lon
     309 
     310    IMPLICIT NONE
     311 
     312      CHARACTER(LEN=*), INTENT(IN) :: field_name
     313      REAL, DIMENSION(:,:,:), INTENT(IN) :: field ! --> field(klon,:,:)
     314 
     315      REAL,DIMENSION(klon_mpi,SIZE(field,2),SIZE(field,3)) :: buffer_omp
     316      REAL :: Field4d(nbp_lon,jj_nb,SIZE(field,2),SIZE(field,3))
     317      INTEGER :: ip, n, nlev, llev
     318 
     319    IF (prt_level >= 10) write(lunout,*)'Begin histrwrite4d_xios ',trim(field_name)
     320 
     321      !Et on.... écrit
     322      IF (SIZE(field,1)/=klon) CALL abort_physic('iophy::histwrite4d','Field first DIMENSION not equal to klon',1)
     323     
     324      nlev=SIZE(field,2)
     325      llev=SIZE(field,3)
     326 
     327      CALL Gather_omp(field,buffer_omp)
     328  !$OMP MASTER
     329      CALL grid1Dto2D_mpi(buffer_omp,field4d)
     330 
     331      CALL xios_send_field(field_name, Field4d(:,:,1:nlev,1:llev))
     332  !$OMP END MASTER   
     333 
     334      IF (prt_level >= 10) write(lunout,*)'End histrwrite4d_xios ',trim(field_name)
     335  END SUBROUTINE histwrite4d_xios
     336 
    293337#endif
    294338
Note: See TracChangeset for help on using the changeset viewer.