Changeset 3090 for trunk/LMDZ.TITAN/libf


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

BdeBatz? : Cleans microphysics and makes few corrections for physics

Location:
trunk/LMDZ.TITAN/libf
Files:
35 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/muphytitan/argparse.F90

    r3083 r3090  
    1 ! Copyright (c) (2013-2015,2017) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
     1! Copyright (c) (2013-2015,2017,2022) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
    22!
    33! This file is part of SWIFT

  • trunk/LMDZ.TITAN/libf/muphytitan/asciiread.f90

    r3083 r3090  
    1 ! Copyright (c) (2013-2015,2017) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
     1! Copyright (c) (2013-2015,2017,2022) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
    22!
    33! This file is part of SWIFT

  • trunk/LMDZ.TITAN/libf/muphytitan/cfgparse.F90

    r3083 r3090  
    1 ! Copyright (c) (2013-2015,2017) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
     1! Copyright (c) (2013-2015,2017,2022) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
    22!
    33! This file is part of SWIFT

  • trunk/LMDZ.TITAN/libf/muphytitan/config/config.ne15.cfg

    r1897 r3090  
    1919#    2 - SF interactions
    2020#    4 - FF interactions.
    21 # (for example: 5 = 4+1 --> SS and FF coagulation only)
     21# (for example : 5 = 4+1 --> SS and FF coagulation only)
    2222haze_coag_interactions = 7
    2323# Enable/disable Haze sedimentation process
    2424haze_sedimentation     = T
    2525# Disable Fiadero correction for sedimentation process
    26 no_fiadero             = F
     26no_fiadero             = T
    2727# Fiadero correction minimum ratio threshold
    2828fiadero_min_ratio      = 0.1
    … …  
    3030fiadero_max_ratio      = 10.
    3131# Force settling velocity to M0
    32 wsed_m0 = T
     32wsed_m0                = T
    3333# Force settling velocity to M3
    34 wsed_m3 = F
     34wsed_m3                = F
    3535# Enable/disable clouds sedimentation process
    3636# (automatically set to F if clouds microphysics is not enabled)
    … …  
    4141# Condensible species configuration file
    4242# (not needed if clouds microphysics is not enabled)
    43 specie_cfg             = ../datagcm/microphysics/mp2m_species.cfg
     43specie_cfg             = ../../datagcm/microphysics/mp2m_species.cfg
    4444
    4545# Enable/disable spherical mode transfert probability
    … …  
    4747# Path of the spherical mode transfert probability look-up tables file
    4848# (optional if 'transfert_probability' is False)
    49 ps2s_file             = ../datagcm/microphysics/mmp_ps2s_rm50_ne15.nc
     49ps2s_file             = ../../datagcm/microphysics/mmp_ps2s_rm50_ne15.nc
    5050
    5151# Electric charging coagulation correction
    … …  
    5454# Path of the electric charging correction factor.
    5555# (optional if 'electric_charging' is False)
    56 mq_file               = ../datagcm/microphysics/mmp_qmean_rm50_ne15.nc
     56mq_file               = ../../datagcm/microphysics/mmp_qmean_rm50_ne15.nc
    5757
    58 #  alpha_X sections contain the parameters of the inter-moments relation function for
    59 #  the mode X: either spherical (s) or fractal (f)
    60 #  dndr_X sections contain the parameters of the size-distribution law of the mode X
     58### Thresholds related parameters
     59### -----------------------------
     60# Aerosol spherical mode total number min. threshold
     61m0as_min = 1e-10
     62# Aerosol spherical mode min. caracteristic radius threshold
     63rcs_min  = 1e-9
     64# Aerosol fractal mode total number min. threshold
     65m0af_min = 1e-10
     66# Aerosol fractal mode min. caracteristic radius threshold (updated to monomer radius at runtime if needed)
     67rcf_min  = 1e-9
     68# cloud drop total number min. threshold
     69m0n_min = 1e-10
     70
     71
     72# =================== #
     73# alpha_k cofficients #
     74# =================== #
     75# alpha_X sections contain the parameters of the inter-moments relation function for the mode X (spherical (s) or fractal (f))
    6176[alpha_s]
    6277a =   1.5044478E-02,  -2.0948806E-01,  -1.5824302E+02,   1.1597818E-01,   9.9502283E-02,  -1.1223796E-01
    6378b =  -2.8622255E-01,   7.7089599E+00,  -1.7000626E+02,   2.6012143E+00,   5.5138784E-01,   9.2024747E-01
    6479c =  -3.0205020E-02,  -3.5510239E+01,  -2.0306468E+02,  -1.3605159E+01,  -4.1653422E+00,  -4.2571698E+00
    65 [dndr_s]
    66 rc = 4.58219580180634588E-007
    67 a0 = 86144.861255561875
    68 c  = 0d0
    69 a  = 2.48333861883769357E-040, 1.46076790655632173E-013, 1.71525517568997062E-009,
    70     1.80855172875974993E-019, 1.48212594918347503E-047, 6.87247318898338451E-081
    71 b  = 59.518212357684796, 15.507500262021228, -5.4179933012448069,
    72     -9.3500794017892854, -18.207927270524777, -27.248924688740562
     80
    7381[alpha_f]
    7482a =   1.5044478E-02,  -2.0948806E-01,  -1.5824302E+02,   1.1597818E-01,   9.9502283E-02,  -1.1223796E-01
    7583b =  -2.8622255E-01,   7.7089599E+00,  -1.7000626E+02,   2.6012143E+00,   5.5138784E-01,   9.2024747E-01
    7684c =  -3.0205020E-02,  -3.5510239E+01,  -2.0306468E+02,  -1.3605159E+01,  -4.1653422E+00,  -4.2571698E+00
    77 [dndr_f]
    78 rc = 4.58219580180634588E-007
    79 a0 = 86144.861255561875
    80 c  = 0d0
    81 a  = 2.48333861883769357E-040, 1.46076790655632173E-013, 1.71525517568997062E-009,
    82     1.80855172875974993E-019, 1.48212594918347503E-047, 6.87247318898338451E-081
    83 b  = 59.518212357684796, 15.507500262021228, -5.4179933012448069,
    84     -9.3500794017892854, -18.207927270524777, -27.248924688740562
    8585
    8686# ================= #
    8787# b^T_k cofficients #
    8888# ================= #
    89 # This section gathers the values of all the btk coefficient used in the coagulation
    90 # equations for the free-molecular regime.
     89# This section gathers the values of all the btk coefficient used in the coagulation equations for the free-molecular regime
    9190[btks]
    9291bt0 = 0.73d0, 0.73d0, 0.75d0, 0.99d0, 0.00d0
    9392bt3 = 0.97d0, 0.97d0, 0.00d0, 0.99d0, 0.99d0
    94 
    95 [optics]
    96 optic_file = /path/to/optics_look_up_table.nc
    97 
    98 
    99 

  • trunk/LMDZ.TITAN/libf/muphytitan/config/mp2m_species.cfg

    r1793 r3090  
    1 #----------------------------------------------------------------------------------------
    2 # mp_scpecies.cfg
     1# ----------------------------------------------------------------------------------------
     2# mp2m_scpecies.cfg
    33# Thermodynamic properties of chemical species used in the cloud microphysics processes.
    44# Each condensible specie that should be treated by the model must be entirely described
    55# here.
    66#
    7 # Air specie is always mandatory. It is saved in a special variable in the model and is
    8 # never used as cloud condensible specie.
    9 #-----------------------------------------------------------------------------------------
     7# ----------------------------------------------------------------------------------------
    108# PARAMETER | DESCRIPTION
    119# ----------|-----------------------------------------------------------------------------
    1210#  name     | Specie name, should be the same as the section name
    13 #  mas      | Molecular weight
    14 #  vol      | Molecular volume
    15 #  ray      | Molecular radius
    16 #  masmol   | Molar mass
    17 #  rho      | Density
    18 #  tc       | Critical temperature
    19 #  tb       | Boiling temperature
    20 #  pc       | Critical pressure (in bar !!!)
     11#  mas      | Molecular weight      [kg]
     12#  vol      | Molecular volume      [m3]
     13#  ray      | Molecular radius      [m]
     14#  masmol   | Molar mass            [kg.mol-1]
     15#  rho      | Density               [kg.m-3]
     16#  tc       | Critical temperature  [K]
     17#  tb       | Boiling temperature   [K]
     18#  pc       | Critical pressure     [bar]
     19#  tx_prod  | Production rate       [kg.m-2.s-1] (actually this is not used by the model)
    2120#  w        | Acentric factor
    22 #  a_sat    | Coefficient A of Psat equation (from Reid et al. 1986).
    23 #  b_sat    | Coefficient B of Psat equation (from Reid et al. 1986).
    24 #  c_sat    | Coefficient C of Psat equation (from Reid et al. 1986).
    25 #  d_sat    | Coefficient D of Psat equation (from Reid et al. 1986).
     21#  a_sat    | Coefficient A of Psat equation (from Reid et al. 1986)
     22#  b_sat    | Coefficient B of Psat equation (from Reid et al. 1986)
     23#  c_sat    | Coefficient C of Psat equation (from Reid et al. 1986)
     24#  d_sat    | Coefficient D of Psat equation (from Reid et al. 1986)
    2625#  mteta    | Wettability (free parameter in most cases, from 0 to 1)
    27 #  tx_prod  | Production rate (actually this is not used by the model)
    28 #-----------------------------------------------------------------------------------------
     26# -----------------------------------------------------------------------------------------
    2927
    30 ### List of actual species to be used in the model:
    31 ### WARNING : the list of species specified here should be ordered:
    32 ###    In the model, ice tracers as well as condensible gazs species must have the same
    33 ###    index.
    34 used_species = "CH4", "C2H6", "C2H2"
    35 
    36 # AIR properties (MANDATORY !!!)
    37 ################################
    38 [air]
    39 name    = "air"
    40 mas     = 4.650e-26
    41 vol     = 5.750e-29
    42 ray     = 1.750e-10
    43 masmol  = 28.e-3
    44 rho     = 808.6
    45 tc      = 126.2
    46 tb      = 77.4
    47 pc      = 33.9
    48 w       = 3.9e-2
    49 a_sat   = -6.09676
    50 b_sat   = 1.13670
    51 c_sat   = -1.04072
    52 d_sat   = -1.93306
    53 mteta   = 0.
    54 tx_prod = 0.
     28### List of actual species to be used in the model :
     29### WARNING : The list of species specified here should be ordered :
     30###           In the model, ice tracers as well as condensible gazs species must
     31###           have the same index.
     32used_species = CH4, C2H6, C2H2, HCN
    5533
    5634# CH4 properties (useful for Titan :)
    5735#####################################
    5836[CH4]
    59 name    = "CH4"
     37name    = CH4
    6038mas     = 2.6578e-26
    6139vol     = 6.252e-29
    6240ray     = 2.000e-10
    63 masmol  = 16.e-3
     41masmol  = 16.04e-3
    6442rho     = 425.
    65 tc      = 190.4
    66 tb      = 111.6
    67 pc      = 46.0
     43tc      = 190.551
     44tb      = 111.66
     45pc      = 45.992
    6846w       = 1.1e-2
    69 a_sat   = -6.00435
    70 b_sat   = 1.18850
    71 c_sat   = -0.83408
    72 d_sat   = -1.22833
    73 mteta   = 0.92
     47a_sat   = -6.02242
     48b_sat   = 1.26652
     49c_sat   = -0.5707
     50d_sat   = -1.366
     51mteta   = 0.99
    7452tx_prod = 0.
    7553
    … …  
    7755#################
    7856[C2H6]
    79 name    = "C2H6"
     57name    = C2H6
    8058mas     = 4.983e-26
    8159vol     = 9.094e-29
    8260ray     = 2.220e-10
    83 masmol  = 30.e-3
     61masmol  = 30.07e-3
    8462rho     = 544.6
    85 tc      = 305.4
    86 tb      = 184.6
    87 pc      = 48.8
     63tc      = 305.33
     64tb      = 184.55
     65pc      = 48.71
    8866w       = 9.9e-2
    89 a_sat   = -6.34307
    90 b_sat   = 1.01163
    91 c_sat   = -1.19116
    92 d_sat   = -2.03539
    93 mteta   = 0.92
     67a_sat   = -6.47500
     68b_sat   = 1.41071
     69c_sat   = -1.1440
     70d_sat   = -1.8590
     71mteta   = 0.99
    9472tx_prod = 1.2e-12
    9573
    … …  
    9775#################
    9876[C2H2]
    99 name    = "C2H2"
     77name    = C2H2
    10078mas     = 4.319e-26
    10179vol     = 7.020e-29
    10280ray     = 2.015e-10
    103 masmol  = 26.e-3
     81masmol  = 26.04e-3
    10482rho     = 615.
    105 tc      = 308.8
    106 tb      = 188.4
    107 pc      = 61.4
     83tc      = 308.35
     84tb      = 188.40
     85pc      = 61.39
    10886w       = 19.0e-2
    109 a_sat   = -6.90128
    110 b_sat   = 1.26873
    111 c_sat   = -2.09113
    112 d_sat   = -2.75601
    113 mteta   = 0.92
     87a_sat   = -6.87886
     88b_sat   = 1.30164
     89c_sat   = -1.22474
     90d_sat   = -3.59556
     91mteta   = 0.99
    11492tx_prod = 3.2e-13
    11593
    … …  
    11795#################
    11896[HCN]
    119 name    = "HCN"
     97name    = HCN
    12098mas     = 4.484e-26
    12199vol     = 6.498e-29
    … …  
    131109c_sat   = -3.004
    132110d_sat   = 1635.
    133 mteta   = 0.92
     111mteta   = 0.99
    134112tx_prod = 1e-12
    135 

  • trunk/LMDZ.TITAN/libf/muphytitan/csystem.c

    r3083 r3090  
    11/*
    2  * Copyright (c) (2013-2015,2017) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
     2 * Copyright (c) (2013-2015,2017,2022) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
    33 *
    44 * This file is part of SWIFT

  • trunk/LMDZ.TITAN/libf/muphytitan/csystem.h

    r3083 r3090  
    11/*
    2  * Copyright (c) (2013-2015,2017) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
     2 * Copyright (c) (2013-2015,2017,2022) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
    33 *
    44 * This file is part of SWIFT

  • trunk/LMDZ.TITAN/libf/muphytitan/defined.h

    r3083 r3090  
    11/*
    2  * Copyright (c) (2013-2015,2017) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
     2 * Copyright (c) (2013-2015,2017,2022) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
    33 *
    44 * This file is part of SWIFT

  • trunk/LMDZ.TITAN/libf/muphytitan/errors.F90

    r3083 r3090  
    1 ! Copyright (c) (2013-2015,2017) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
     1! Copyright (c) (2013-2015,2017,2022) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
    22!
    33! This file is part of SWIFT

  • trunk/LMDZ.TITAN/libf/muphytitan/fsystem.F90

    r3083 r3090  
    1 ! Copyright (c) (2013-2015,2017) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
     1! Copyright (c) (2013-2015,2017,2022) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
    22!
    33! This file is part of SWIFT

  • trunk/LMDZ.TITAN/libf/muphytitan/mm_clouds.f90

    r3083 r3090  
    1 ! Copyright 2013-2015,2017,2022 Université de Reims Champagne-Ardenne
    2 ! Contributors: J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
     1! Copyright (2013-2015,2017,2022-2023) Université de Reims Champagne-Ardenne
     2! Contributors : J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
    33! email of the author : jeremie.burgalat@univ-reims.fr
    44!
    … …  
    3535!! summary: Clouds microphysics module
    3636!! author: J. Burgalat
    37 !! date: 2013-2015,2017,2022
    38 !! corrections: B. de Batz de Trenquelléon (2023)
     37!! date: 2013-2015,2017,2022-2023
     38!! modifications: B. de Batz de Trenquelléon
    3939
    4040MODULE MM_CLOUDS
    … …  
    109109      call mm_cloud_sedimentation(zdm0n,zdm3n,zdm3i)
    110110
    111       ! computes precipitation, settling velocity and flux of ices
     111      ! Computes settling velocity [m.s-1] of clouds (ccn and ices)
     112      mm_ccn_vsed(:) = wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad)
     113
     114      ! Computes flux [kg.m-2.s-1] and precipitation [m.iphysiq] of ccn
     115      mm_ccn_flux(:) = get_mass_flux(mm_rhoaer,mm_m3ccn(:))
    112116      mm_ccn_prec = SUM(zdm3n*mm_dzlev)
    113       mm_ccn_w(:) = wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad)
    114       mm_ccn_flux(:) = get_mass_flux(mm_rhoaer,mm_m3ccn(:))
    115 
    116       DO i=1, mm_nesp
     117
     118      ! Computes flux [kg.m-2.s-1] and precipitation [m.iphysiq] of ices
     119      DO i = 1, mm_nesp
     120        mm_ice_fluxes(:,i) = get_mass_flux(mm_xESPS(i)%rho,(3._mm_wp*mm_m3ice(:,i))/(4._mm_wp*mm_pi))
    117121        mm_ice_prec(i) = SUM(zdm3i(:,i)*mm_dzlev)
    118         mm_ice_fluxes(:,i) = get_mass_flux(mm_xESPS(i)%rho,mm_m3ice(:,i))
    119122      ENDDO
     123
    120124      ! updates tendencies
    121125      dm0n = dm0n + zdm0n
    … …  
    724728    REAL(kind=mm_wp), INTENT(in) :: rad !! Radius of the particle (m).
    725729    REAL(kind=mm_wp) :: w               !! Settling velocity (\(m.s^{-1}\)).
    726     REAL(kind=mm_wp) :: Us, Fc
     730    REAL(kind=mm_wp) :: Us, Fc, kn
     731    REAL(kind=mm_wp), PARAMETER :: ra = 1.75e-10_mm_wp
     732   
     733    ! Knudsen number
     734    kn = (mm_kboltz * t) / (p * 4._mm_wp * sqrt(2._mm_wp) * mm_pi * ra**2) / rad
     735
     736    ! Computes slip-flow correction : Fc = 1 + (mm_akn * mm_lambda_g(t,p) / rad)
     737    Fc = (1._mm_wp + 1.2517_mm_wp*kn + 0.4_mm_wp*kn*dexp(-1.1_mm_wp/kn))
    727738   
    728739    ! Computes Stokes settling velocity
    729740    Us = (2._mm_wp * rad**2 * rho * mm_effg(z)) / (9._mm_wp * mm_eta_g(t))
    730741   
    731     ! Apply slip-flow correction
    732     Fc = 1 + (mm_akn * mm_lambda_g(t,p) / rad)
    733 
    734     ! Computes settling velocity
    735     w = Us * Fc
    736     !>>> [TEMPO : BBT]
    737     !w = Us * Fc * 5.
    738     !<<< [TEMPO : BBT]
     742    ! Computes settling velocity (correction factor : x2.0)
     743    w = Us * Fc * 2._mm_wp
    739744  END FUNCTION wsettle
    740745
    … …  
    753758    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: m3
    754759    !! Vertical profile of the total volume of tracer (i.e. M3) from __TOP__ to __GROUND__ (\(m^{3}.m^{-3}\)).
    755     REAL(kind=mm_wp), DIMENSION(SIZE(m3)) :: flx
     760    REAL(kind=mm_wp), DIMENSION(SIZE(m3))      :: flx
    756761    !! Mass sedimentation fluxes at each layer from __TOP__ to __GROUND__ (\(kg.m^{-2}.s^{-1}\)).
    757     REAL(kind=mm_wp), SAVE :: fac = 4._mm_wp/3._mm_wp * mm_pi
    758     flx = fac * rho * m3 * wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad)
     762    REAL(kind=mm_wp), SAVE :: pifac = (4._mm_wp * mm_pi) / 3._mm_wp
     763    flx = rho * pifac * m3 * wsettle(mm_play,mm_temp,mm_zlay,mm_drho,mm_drad)
    759764    RETURN
    760765  END FUNCTION get_mass_flux

  • trunk/LMDZ.TITAN/libf/muphytitan/mm_globals.f90

    r3083 r3090  
    1 ! Copyright 2013-2015,2017,2022 Université de Reims Champagne-Ardenne
    2 ! Contributor: J. Burgalat (GSMA, URCA)
     1! Copyright (2013-2015,2017,2022-2023) Université de Reims Champagne-Ardenne
     2! Contributors : J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
    33! email of the author : jeremie.burgalat@univ-reims.fr
    44!
    … …  
    3535!! summary: Parameters and global variables module.
    3636!! author: J. Burgalat
    37 !! date: 2013-2015,2017,2022
     37!! date: 2013-2015,2017,2022-2023
     38!! modifications: B. de Batz de Trenquelléon
    3839
    3940MODULE MM_GLOBALS
    … …  
    274275
    275276  ! Thresholds !
    276 
    277277  !> (min.) Total number of aerosols minimum threshold for the spherical mode.
    278278  REAL(kind=mm_wp), SAVE :: mm_m0as_min = 1.e-10_mm_wp
    279 
    280279  !> (min.) Total volume of aerosols minimum threshold for the spherical mode.
    281280  REAL(kind=mm_wp), SAVE :: mm_m3as_min = 1.e-40_mm_wp
    282 
    283281  !> Characteristic radius minimum threshold for the spherical mode.
    284282  REAL(kind=mm_wp), SAVE :: mm_rcs_min = 1.e-9_mm_wp
    … …  
    286284  !> (min.) Total number of aerosols minimum threshold for the fractal mode.
    287285  REAL(kind=mm_wp), SAVE :: mm_m0af_min = 1.e-10_mm_wp
    288 
    289286  !> (min.) Total volume of aerosols minimum threshold for the fractal mode.
    290287  REAL(kind=mm_wp), SAVE :: mm_m3af_min = 1.e-40_mm_wp
    291 
    292288  !> Characteristic radius minimum threshold for the fractal mode.
    293289  REAL(kind=mm_wp), SAVE :: mm_rcf_min = 1.e-9_mm_wp
    … …  
    295291  !> (min.) Total number of cloud drop minimum threshold.
    296292  REAL(kind=mm_wp), SAVE :: mm_m0n_min = 1.e-10_mm_wp
    297 
    298293  !> (min.) Total volume of cloud drop minimum threshold.
    299294  REAL(kind=mm_wp), SAVE :: mm_m3cld_min = 1.e-40_mm_wp
     295  !> Minimum cloud drop radius
     296  REAL(kind=mm_wp), SAVE :: mm_drad_min = 1.e-9_mm_wp
     297  !> Maximum cloud drop radius
     298  REAL(kind=mm_wp), SAVE :: mm_drad_max = 1.e-3_mm_wp
    300299
    301300  !> Characteristic radius threshold.
    302301  REAL(kind=mm_wp), SAVE :: mm_rc_min = 1.e-200_mm_wp
    303 
    304   !> Minimum cloud drop radius
    305   REAL(kind=mm_wp), SAVE :: mm_drad_min = 1.e-9_mm_wp
    306 
    307   !> Maximum cloud drop radius
    308   REAL(kind=mm_wp), SAVE :: mm_drad_max = 1.e-3_mm_wp
    309302
    310303  !> Name of condensible species.
    … …  
    442435  REAL(kind=mm_wp), SAVE :: mm_aer_prec = 0._mm_wp
    443436
     437  !> CCN precipitations (m).
     438  !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
     439  REAL(kind=mm_wp), SAVE :: mm_ccn_prec = 0._mm_wp
     440
    444441  !> Spherical mode \(M_{0}\) settling velocity (\(m.s^{-1}\)).
    445442  !!
    … …  
    478475  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_m3af_vsed
    479476
     477  !> CCN settling velocity (\(m.s^{-1}\)).
     478  !!
     479  !! It is a vector with the vertical layers that contains the
     480  !! settling velocity for CCN (and ices).
     481  !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
     482  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ccn_vsed
     483
    480484  !> Spherical aerosol mass fluxes (\(kg.m^{-2}.s^{-1}\)).
    481485  !!
    … …  
    493497  !! This variable is always negative.
    494498  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_aer_f_flux
    495 
    496   !> CCN precipitations (m).
    497   !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
    498   REAL(kind=mm_wp), SAVE :: mm_ccn_prec = 0._mm_wp
    499 
    500   !> CCN settling velocity (\(m.s^{-1}\)).
    501   !!
    502   !! It is a vector with the vertical layers that contains the
    503   !! settling velocity for CCN (and ices).
    504   !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
    505   !! @note
    506   !! This variable is always positive.
    507   REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ccn_w
    508499 
    509500  !> CCN mass fluxes (\(kg.m^{-2}.s^{-1}\)).
    … …  
    512503  !! mass fluxes for CCN.
    513504  !! It is updated in [[mm_clouds(module):mm_cloud_microphysics(subroutine)]].
    514   !! @note
    515   !! This variable is always positive.
    516505  REAL(kind=mm_wp), DIMENSION(:), ALLOCATABLE, SAVE :: mm_ccn_flux
    517506
    … …  
    593582  !$OMP THREADPRIVATE(mm_m0ccn,mm_m3ccn,mm_m3ice,mm_gazs)
    594583  !$OMP THREADPRIVATE(mm_rcs,mm_rcf,mm_drad,mm_drho)
    595   !$OMP THREADPRIVATE(mm_aer_s_flux,mm_aer_f_flux,mm_ccn_w,mm_ccn_flux,mm_ice_prec,mm_ice_fluxes,mm_gazs_sat)
    596584  !$OMP THREADPRIVATE(mm_m0as_vsed,mm_m3as_vsed,mm_m0af_vsed,mm_m3af_vsed)
     585  !$OMP THREADPRIVATE(mm_aer_s_flux,mm_aer_f_flux,mm_ccn_vsed,mm_ccn_flux,mm_ice_prec,mm_ice_fluxes,mm_gazs_sat)
    597586  !$OMP THREADPRIVATE(mm_m0as_min,mm_m3as_min,mm_rcs_min,mm_m0af_min,mm_m3af_min,mm_rcf_min,mm_m0n_min,mm_m3cld_min)
    598587  !$OMP THREADPRIVATE(mm_nla,mm_nle)
    … …  
    997986    err = mm_check_opt(cfg_get_value(cfg,"fiadero_min_ratio",mm_fiadero_min),mm_fiadero_min,zfiamin,wlog=mm_log)
    998987    err = mm_check_opt(cfg_get_value(cfg,"fiadero_max_ratio",mm_fiadero_max),mm_fiadero_max,zfiamax,wlog=mm_log)
    999 
    1000988    err = mm_check_opt(cfg_get_value(cfg,"m0as_min",mm_m0as_min),mm_m0as_min,1e-10_mm_wp,wlog=mm_log)
    1001989    err = mm_check_opt(cfg_get_value(cfg,"rcs_min",mm_rcs_min),mm_rcs_min,1e-9_mm_wp,wlog=mm_log)
    … …  
    1003991    err = mm_check_opt(cfg_get_value(cfg,"rcf_min",mm_rcf_min),mm_rcf_min,mm_rm,wlog=mm_log)
    1004992    err = mm_check_opt(cfg_get_value(cfg,"m0n_min",mm_m0n_min),mm_m0n_min,1e-10_mm_wp,wlog=mm_log)
    1005 
    1006993
    1007994    ! force fractal radius minimum threshold to monomer radius ^^
    … …  
    11521139    IF (.NOT.ALLOCATED(mm_rcf))     ALLOCATE(mm_rcf(mm_nla))
    11531140    ! Allocate memory for diagnostics
     1141    IF (.NOT.ALLOCATED(mm_m0as_vsed)) THEN
     1142      ALLOCATE(mm_m0as_vsed(mm_nla)) ; mm_m0as_vsed(:) = 0._mm_wp
     1143    ENDIF
     1144    IF (.NOT.ALLOCATED(mm_m3as_vsed)) THEN
     1145      ALLOCATE(mm_m3as_vsed(mm_nla)) ; mm_m3as_vsed(:) = 0._mm_wp
     1146    ENDIF
     1147    IF (.NOT.ALLOCATED(mm_m0af_vsed)) THEN
     1148      ALLOCATE(mm_m0af_vsed(mm_nla)) ; mm_m0af_vsed(:) = 0._mm_wp
     1149    ENDIF
     1150    IF (.NOT.ALLOCATED(mm_m3af_vsed)) THEN
     1151      ALLOCATE(mm_m3af_vsed(mm_nla)) ; mm_m3af_vsed(:) = 0._mm_wp
     1152    ENDIF
    11541153    IF (.NOT.ALLOCATED(mm_aer_s_flux)) THEN
    11551154      ALLOCATE(mm_aer_s_flux(mm_nla)) ; mm_aer_s_flux(:) = 0._mm_wp
    … …  
    11571156    IF (.NOT.ALLOCATED(mm_aer_f_flux)) THEN
    11581157      ALLOCATE(mm_aer_f_flux(mm_nla)) ; mm_aer_f_flux(:) = 0._mm_wp
    1159     ENDIF
    1160     IF (.NOT.ALLOCATED(mm_m0as_vsed)) THEN
    1161       ALLOCATE(mm_m0as_vsed(mm_nla)) ; mm_m0as_vsed(:) = 0._mm_wp
    1162     ENDIF
    1163     IF (.NOT.ALLOCATED(mm_m3as_vsed)) THEN
    1164       ALLOCATE(mm_m3as_vsed(mm_nla)) ; mm_m3as_vsed(:) = 0._mm_wp
    1165     ENDIF
    1166     IF (.NOT.ALLOCATED(mm_m0af_vsed)) THEN
    1167       ALLOCATE(mm_m0af_vsed(mm_nla)) ; mm_m0af_vsed(:) = 0._mm_wp
    1168     ENDIF
    1169     IF (.NOT.ALLOCATED(mm_m3af_vsed)) THEN
    1170       ALLOCATE(mm_m3af_vsed(mm_nla)) ; mm_m3af_vsed(:) = 0._mm_wp
    11711158    ENDIF
    11721159    ! note : mm_dzlev is already from top to ground
    … …  
    12411228    IF (.NOT.ALLOCATED(mm_drho))    ALLOCATE(mm_drho(mm_nla))
    12421229    ! Allocate memory for diagnostics
    1243     IF (.NOT.ALLOCATED(mm_ccn_w)) THEN
    1244       ALLOCATE(mm_ccn_w(mm_nla)) ; mm_ccn_w(:) = 0._mm_wp
     1230    IF (.NOT.ALLOCATED(mm_ccn_vsed)) THEN
     1231      ALLOCATE(mm_ccn_vsed(mm_nla)) ; mm_ccn_vsed(:) = 0._mm_wp
    12451232    ENDIF
    12461233    IF (.NOT.ALLOCATED(mm_ccn_flux)) THEN
    … …  
    14441431    ! Initialization :
    14451432    Ntot = m0ccn
    1446     Vtot = pifac * m3ccn + SUM(m3ice)
    1447     Wtot = pifac * ((m3ccn*mm_rhoaer) + SUM(m3ice*mm_xESPS(:)%rho))
     1433    Vtot = pifac*m3ccn + SUM(m3ice)
     1434    Wtot = pifac*m3ccn*mm_rhoaer + SUM(m3ice*mm_xESPS(:)%rho)
    14481435
    14491436    IF (Ntot <= mm_m0n_min .OR. Vtot <= mm_m3cld_min) THEN

  • trunk/LMDZ.TITAN/libf/muphytitan/mm_haze.f90

    r3083 r3090  
    731731      cslf = mm_akn * mm_lambda_g(mm_btemp(i),mm_plev(i))
    732732      wth(i) = - csto/(rb2ra*afun(k)) * (rc(i-1)**ar1 * afun(af1) + cslf/rb2ra * rc(i-1)**ar2 * afun(af2))
    733 
    734       ! >>> [TEMPO : BBT]
    735       !wth(i) = wth(i) * (2574e3 / (2574e3+mm_zlev(i)))**4
    736       ! <<< [TEMPO : BBT]
    737733     
    738734      ! now correct velocity to reduce numerical diffusion

  • trunk/LMDZ.TITAN/libf/muphytitan/mm_methods.f90

    r3083 r3090  
    1 ! Copyright 2013-2015,2017,2022 Université de Reims Champagne-Ardenne
    2 ! Contributors: J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
     1! Copyright (2013-2015,2017,2022-2023) Université de Reims Champagne-Ardenne
     2! Contributors : J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
    33! email of the author : jeremie.burgalat@univ-reims.fr
    44!
    … …  
    3535!! summary: Model miscellaneous methods module.
    3636!! author: J. Burgalat
    37 !! date: 2013-2015,2017,2022
    38 !! corrections: B. de Batz de Trenquelléon (2023)
     37!! date: 2013-2015,2017,2022-2023
     38!! modifications: B. de Batz de Trenquelléon
    3939
    4040MODULE MM_METHODS
    … …  
    133133  !!
    134134  !! The method computes the molar mixing ratio at saturation of a given specie at given temperature(s)
    135   !! and pressure level(s) [Fray and Schmidt (2009)].
     135  !! and pressure level(s).
    136136  !!
    137137  !! ```fortran
    … …  
    314314    !!
    315315    !! @warning
    316     !! The method applies a multiplicative factor of 0.85 if the specie is CH4 :
     316    !! The method applies a multiplicative factor of 0.80 if the specie is CH4 :
    317317    !! this is done to account for dissolution in N2 and is somehow specific to Titan atmosphere.
    318318    REAL(kind=mm_wp), INTENT(in) :: temp !! Temperature (K).
    319319    REAL(kind=mm_wp), INTENT(in) :: pres !! Pressure level (Pa).
    320     TYPE(mm_esp), INTENT(in)    :: xESP !! Specie properties.
    321     REAL(kind=mm_wp) :: res              !! Mass mixing ratio of the specie.
    322     REAL(kind=mm_wp) :: psat
     320    TYPE(mm_esp), INTENT(in)     :: xESP !! Specie properties.
     321    REAL(kind=mm_wp)             :: res  !! Mass mixing ratio of the specie.
     322    REAL(kind=mm_wp)             :: psat
    323323    psat = mm_psatX(temp,xESP)
    324324    res = (psat / pres) * xESP%fmol2fmas
    325     ! Peculiar case : CH4 : x0.85 (dissolution in N2)
     325    ! Peculiar case of CH4 : x0.80 (dissolution in N2)
    326326    IF (xESP%name == "CH4") THEN
    327       res = res * 0.85_mm_wp
     327      res = res * 0.80_mm_wp
    328328      IF (mm_debug) WRITE(*,'(a)') "[DEBUG] mm_qsat: applying .85 factor to qsat for CH4 specie (N2 dissolution)"
    329329    ENDIF
    … …  
    339339    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: pres !! Pressure levels (Pa).
    340340    TYPE(mm_esp), INTENT(in)                   :: xESP !! Specie properties.
    341     REAL(kind=mm_wp), DIMENSION(SIZE(temp)) :: res     !! Mass mixing ratios of the specie.
    342     REAL(kind=mm_wp), DIMENSION(SIZE(temp)) :: psat
     341    REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  !! Mass mixing ratios of the specie.
     342    REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: psat
    343343    psat = mm_psatX(temp,xESP)
    344344    res = (psat / pres) * xESP%fmol2fmas
    345     ! Peculiar case : CH4 : x0.85 (dissolution in N2)
     345    ! Peculiar case of CH4 : x0.80 (dissolution in N2)
    346346    IF (xESP%name == "CH4") THEN
    347       res = res * 0.85_mm_wp
     347      res = res * 0.80_mm_wp
    348348      IF (mm_debug) WRITE(*,'(a)') "[DEBUG] mm_qsat: applying .85 factor to qsat for CH4 specie (N2 dissolution)"
    349349    ENDIF
    … …  
    352352  FUNCTION ysatX_sc(temp,pres,xESP) RESULT(res)
    353353    !! Get the molar mixing ratio of a given specie at saturation (scalar).
     354    !! Compute saturation profiles (mol/mol) for condensable tracers
    354355    !!
    355356    !! @warning
    356     !! The method applies a multiplicative factor of 0.85 if the specie is CH4 :
     357    !! The method applies a multiplicative factor of 0.80 if the specie is CH4 :
    357358    !! this is done to account for dissolution in N2 and is somehow specific to Titan atmosphere.
    358359    REAL(kind=mm_wp), INTENT(in) :: temp !! Temperature (K).
    … …  
    361362    REAL(kind=mm_wp)             :: res  !! Molar mixing ratio of the specie.
    362363
    363     ! Fray and Schmidt (2009)
    364364    IF(xESP%name == "C2H2") THEN
     365      ! Fray and Schmidt (2009)
    365366      res = (1.0e5 / pres) * exp(1.340e1 - 2.536e3/temp)
    366367   
    367368    ELSE IF(xESP%name == "C2H6") THEN
     369      ! Fray and Schmidt (2009)
    368370      res = (1.0e5 / pres) * exp(1.511e1 - 2.207e3/temp - 2.411e4/temp**2 + 7.744e5/temp**3 - 1.161e7/temp**4 + 6.763e7/temp**5)
    369371   
    370372    ELSE IF(xESP%name == "HCN") THEN
     373      ! Fray and Schmidt (2009)
    371374      res = (1.0e5 / pres) * exp(1.393e1 - 3.624e3/temp - 1.325e5/temp**2 + 6.314e6/temp**3 - 1.128e8/temp**4)
    372375   
    373376    ELSE IF (xESP%name == "CH4") THEN
     377      ! Fray and Schmidt (2009)
    374378      res = (1.0e5 / pres) * exp(1.051e1 - 1.110e3/temp - 4.341e3/temp**2 + 1.035e5/temp**3 - 7.910e5/temp**4)
    375       res = res * 0.85_mm_wp
    376       !IF (res < 0.014) THEN
    377       !  res = 0.014
    378       !ENDIF
     379      ! Peculiar case of CH4 : x0.80 (dissolution in N2)
     380      res = res * 0.80_mm_wp
     381      ! Forcing CH4 to 1.4% minimum
     382      IF (res < 0.014) THEN
     383        res = 0.014
     384      ENDIF
    379385    ENDIF
    380386  END FUNCTION ysatX_sc
    … …  
    382388  FUNCTION ysatX_ve(temp,pres,xESP) RESULT(res)
    383389    !! Get the molar mixing ratio of a given specie at saturation (vector).
     390    !! Compute saturation profiles (mol/mol) for condensable tracers
    384391    !!
    385392    !! @warning
    386     !! The method applies a multiplicative factor of 0.85 if the specie is CH4 :
     393    !! The method applies a multiplicative factor of 0.80 if the specie is CH4 :
    387394    !! this is done to account for dissolution in N2 and is somehow specific to Titan atmosphere.
    388395    REAL(kind=mm_wp), INTENT(in), DIMENSION(:) :: temp !! Temperatures (K).
    … …  
    391398    REAL(kind=mm_wp), DIMENSION(SIZE(temp))    :: res  !! Molar mixing ratios of the specie.
    392399
    393     ! Fray and Schmidt (2009)
    394400    IF(xESP%name == "C2H2") THEN
     401      ! Fray and Schmidt (2009)
    395402      res = (1.0e5 / pres) * exp(1.340e1 - 2.536e3/temp)
    396403   
    397404    ELSE IF(xESP%name == "C2H6") THEN
     405      ! Fray and Schmidt (2009)
    398406      res = (1.0e5 / pres) * exp(1.511e1 - 2.207e3/temp - 2.411e4/temp**2 + 7.744e5/temp**3 - 1.161e7/temp**4 + 6.763e7/temp**5)
    399407   
    400408    ELSE IF(xESP%name == "HCN") THEN
     409      ! Fray and Schmidt (2009)
    401410      res = (1.0e5 / pres) * exp(1.393e1 - 3.624e3/temp - 1.325e5/temp**2 + 6.314e6/temp**3 - 1.128e8/temp**4)
    402411   
    … …  
    404413    ELSE IF (xESP%name == "CH4") THEN
    405414      res = (1.0e5 / pres) * exp(1.051e1 - 1.110e3/temp - 4.341e3/temp**2 + 1.035e5/temp**3 - 7.910e5/temp**4)
    406       res = res * 0.85_mm_wp
    407       !WHERE (res(:) < 0.014) res(:) = 0.014
     415      ! Peculiar case of CH4 : x0.80 (dissolution in N2)
     416      res = res * 0.80_mm_wp
     417      ! Forcing CH4 to 1.4% minimum
     418      WHERE (res(:) < 0.014) res(:) = 0.014
    408419    ENDIF
    409420  END FUNCTION ysatX_ve

  • trunk/LMDZ.TITAN/libf/muphytitan/mm_microphysic.f90

    r3083 r3090  
    1 ! Copyright 2013-2015,2017,2022 Université de Reims Champagne-Ardenne
    2 ! Contributor: J. Burgalat (GSMA, URCA)
     1! Copyright (2013-2015,2017,2022-2023) Université de Reims Champagne-Ardenne
     2! Contributors : J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
    33! email of the author : jeremie.burgalat@univ-reims.fr
    44!
    … …  
    3535!! brief: Microphysic processes interface module.
    3636!! author: J. Burgalat
    37 !! date: 2013-2015,2017,2022
     37!! date: 2013-2015,2017,2022-2023
     38!! modifications: B. de Batz de Trenquelléon
    3839
    3940MODULE MM_MICROPHYSIC
    … …  
    114115      dm0n   = dm0n(mm_nla:1:-1) * mm_dzlev(mm_nla:1:-1)
    115116      dm3n   = dm3n(mm_nla:1:-1) * mm_dzlev(mm_nla:1:-1)
    116       ! sanity check for clouds tendencies : call mm_check_tendencies(mm_m0ccn,dm0n)
    117       ! sanity check for clouds tendencies : call mm_check_tendencies(mm_m3ccn,dm3n)
    118117      DO i=1,mm_nesp
    119118        dm3i(:,i)  = dm3i(mm_nla:1:-1,i) * mm_dzlev(mm_nla:1:-1)
    120         ! sanity check for clouds tendencies : call mm_check_tendencies(mm_m3ice,dm3i)
    121119        dgazs(:,i) = dgazs(mm_nla:1:-1,i)
    122         ! no sanity check for gazs, let's prey.
    123120      ENDDO
    124121    ELSE
    … …  
    130127    dm0a_f = dm0a_f(mm_nla:1:-1) * mm_dzlev(mm_nla:1:-1)
    131128    dm3a_f = dm3a_f(mm_nla:1:-1) * mm_dzlev(mm_nla:1:-1)
    132     ! sanity check
    133     ! sanity check : call mm_check_tendencies(mm_m0aer_s,dm0a_s)
    134     ! sanity check : call mm_check_tendencies(mm_m3aer_s,dm3a_s)
    135     ! sanity check : call mm_check_tendencies(mm_m0aer_f,dm0a_f)
    136     ! sanity check : call mm_check_tendencies(mm_m3aer_f,dm3a_f)
    137129    RETURN
    138130  END FUNCTION muphys_all
    … …  
    168160    dm0a_f = dm0a_f(mm_nla:1:-1) * mm_dzlev(mm_nla:1:-1)
    169161    dm3a_f = dm3a_f(mm_nla:1:-1) * mm_dzlev(mm_nla:1:-1)
    170     ! sanity check
    171     ! sanity check : call mm_check_tendencies(mm_m0aer_s,dm0a_s)
    172     ! sanity check : call mm_check_tendencies(mm_m3aer_s,dm3a_s)
    173     ! sanity check : call mm_check_tendencies(mm_m0aer_f,dm0a_f)
    174     ! sanity check : call mm_check_tendencies(mm_m3aer_f,dm3a_f)
    175162    RETURN
    176163  END FUNCTION muphys_nocld
    … …  
    219206      IF (PRESENT(ccn_prec))   ccn_prec   = ABS(mm_ccn_prec)
    220207      IF (PRESENT(ice_prec))   ice_prec   = ABS(mm_ice_prec)
    221       IF (PRESENT(ccn_w))      ccn_w      = mm_ccn_w(mm_nla:1:-1)
    222       IF (PRESENT(ccn_flux))   ccn_flux   = -mm_ccn_flux(mm_nla:1:-1)
     208      IF (PRESENT(ccn_w))      ccn_w      = mm_ccn_vsed(mm_nla:1:-1)
     209      IF (PRESENT(ccn_flux))   ccn_flux   = mm_ccn_flux(mm_nla:1:-1)
    223210      IF (PRESENT(ice_fluxes)) ice_fluxes = mm_ice_fluxes(mm_nla:1:-1,:)
    224       IF (PRESENT(gazs_sat))   gazs_sat   =  mm_gazs_sat(mm_nla:1:-1,:)
     211      IF (PRESENT(gazs_sat))   gazs_sat   = mm_gazs_sat(mm_nla:1:-1,:)
    225212    ELSE
    226213      IF (PRESENT(ccn_prec))   ccn_prec   = 0._mm_wp

  • trunk/LMDZ.TITAN/libf/muphytitan/mmp_gcm.f90

    r3083 r3090  
    1 ! Copyright 2017 Université de Reims Champagne-Ardenne
    2 ! Contributor: J. Burgalat (GSMA, URCA)
     1! Copyright (2017,2022-2023) Université de Reims Champagne-Ardenne
     2! Contributors : J. Burgalat (GSMA, URCA), B. de Batz de Trenquelléon (GSMA, URCA)
    33! email of the author : jeremie.burgalat@univ-reims.fr
    44!
    … …  
    3535!! summary: YAMMS interfaces for the LMDZ GCM.
    3636!! author: J. Burgalat
    37 !! date: 2017
     37!! date: 2017,2022,2023
     38!! modifications: B. de Batz de Trenquelléon
    3839MODULE MMP_GCM
    3940  !! Interface to YAMMS for the LMDZ GCM.
    … …  
    247248    WRITE(*,'(a,L2)')     "electric_charging    : ", mmp_w_qe
    248249    call mm_dump_parameters()
    249 
    250250    IF (clouds) THEN
    251251      DO i=1, size(mm_xESPS)

  • trunk/LMDZ.TITAN/libf/muphytitan/string_op.F90

    r3083 r3090  
    1 ! Copyright (c) (2013-2015,2017) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
     1! Copyright (c) (2013-2015,2017,2022) Jeremie Burgalat (jeremie.burgalat@univ-reims.fr).
    22!
    33! This file is part of SWIFT

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