Changeset 4601 for LMDZ6


Ignore:
Timestamp:
Jul 1, 2023, 12:07:30 AM (17 months ago)
Author:
dcugnet
Message:

StratAer? commit, N. Lebas

Location:
LMDZ6/trunk
Files:
5 added
1 deleted
15 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/DefLists/field_def_lmdz.xml

    r4575 r4601  
    749749        <field id="heat_volc"    long_name="SW heating rate from volcano"    unit="K/s" />
    750750        <field id="cool_volc"    long_name="LW cooling rate from volcano"    unit="K/s" />
     751        <field id="dqch4"        long_name="H2O due to CH4 oxidation and photolysis" unit="(kg/kg)/s" />
    751752        <field id="pvap"    long_name="pvap intermediary variable" unit="-">pres*ovap*461.5 / (287.04*(1.+ (10.9491/18.0153)*ovap)) </field>
    752753        <field id="oclr"    long_name="Clear sky total water"    unit="kg/kg" />

  • LMDZ6/trunk/DefLists/file_def_histdaystrataer_lmdz.xml

    r3533 r4601  
    11<file_definition>
    2     <file_group id="defile">
    3         <file id="histdaystrataer" name="histdaystrataer" output_freq="1d" output_level="5" enabled="_AUTO_" compression_level="4">
     2  <file_group id="defile">
     3    <file id="histdaystrataer" name="histdaystrataer" output_freq="1d" output_level="_AUTO_" enabled="_AUTO_" compression_level="4">
     4     
     5      <field_group grid_ref="grid_out"  level="3">
     6        <field field_ref="OD550_strat_only" level="1" />
     7        <field field_ref="OD1020_strat_only" level="1" />
     8        <field field_ref="surf_PM25_sulf" level="2" />
     9        <field field_ref="budg_dep_dry_ocs" level="3" />
     10        <field field_ref="budg_dep_wet_ocs" level="3" />
     11        <field field_ref="budg_dep_dry_so2" level="3" />
     12        <field field_ref="budg_dep_wet_so2" level="3" />
     13        <field field_ref="budg_dep_dry_h2so4" level="3" />
     14        <field field_ref="budg_dep_wet_h2so4" level="3" />
     15        <field field_ref="budg_dep_dry_part" level="3" />
     16        <field field_ref="budg_dep_wet_part" level="3" />
     17        <field field_ref="budg_emi_ocs" level="3" />
     18        <field field_ref="budg_emi_so2" level="3" />
     19        <field field_ref="budg_emi_h2so4" level="3" />
     20        <field field_ref="budg_emi_part" level="3" />
     21        <field field_ref="budg_ocs_to_so2" level="3" />
     22        <field field_ref="budg_so2_to_h2so4" level="3" />
     23        <field field_ref="budg_h2so4_to_part" level="3" />
     24        <field field_ref="budg_sed_part" level="3" />
     25        <field field_ref="p_tropopause" level="1" />
     26        <field field_ref="z_tropopause" level="1" />
     27        <field field_ref="t_tropopause" level="3" />
     28      </field_group>
     29     
     30      <field_group grid_ref="grid_out_presnivs" level="10">
     31        <field field_ref="ext_strat_550" level="1" />
     32        <field field_ref="ext_strat_1020" level="5" />
     33        <field field_ref="budg_3D_nucl" level="10" />
     34        <field field_ref="budg_3D_cond_evap" level="10" />
     35        <field field_ref="budg_3D_ocs_to_so2" level="10" />
     36        <field field_ref="budg_3D_so2_to_h2so4" level="10" />
     37        <field field_ref="budg_3D_backgr_ocs"  level="10" />
     38        <field field_ref="budg_3D_backgr_so2" level="10" />
     39        <field field_ref="R2SO4" level="1" />
     40        <field field_ref="OCS_lifetime" level="10" />
     41        <field field_ref="SO2_lifetime" level="10" />
     42        <field field_ref="vsed_aer" level="10" />
     43        <field field_ref="f_r_wet" level="2" />
     44        <field field_ref="mass"  level="10" />
     45        <field field_ref="temp"  level="10" />
     46        <field field_ref="pres"    level="10" />
     47        <field field_ref="h2o"     level="1" />
     48      </field_group>
     49     
     50      <field_group grid_ref="grid_out_presnivs" level="10">
     51        <field field_ref="GASOCS"  level="1" />
     52        <field field_ref="GASSO2"  level="1" />
     53        <field field_ref="GASH2SO4"  level="1" />
     54        <field field_ref="BIN01" level="2" />
     55        <field field_ref="BIN02" level="2" />
     56        <field field_ref="BIN03" level="2" />
     57        <field field_ref="BIN04" level="2" />
     58        <field field_ref="BIN05" level="2" />
     59        <field field_ref="BIN06" level="2" />
     60        <field field_ref="BIN07" level="2" />
     61        <field field_ref="BIN08" level="2" />
     62        <field field_ref="BIN09" level="2" />
     63        <field field_ref="BIN10" level="2" />
     64        <field field_ref="BIN11" level="2" />
     65        <field field_ref="BIN12" level="2" />
     66        <field field_ref="BIN13" level="2" />
     67        <field field_ref="BIN14" level="2" />
     68        <field field_ref="BIN15" level="2" />
     69        <field field_ref="BIN16" level="2" />
     70        <field field_ref="BIN17" level="2" />
     71        <field field_ref="BIN18" level="2" />
     72        <field field_ref="BIN19" level="2" />
     73        <field field_ref="BIN20" level="2" />
     74        <field field_ref="BIN21" level="2" />
     75        <field field_ref="BIN22" level="2" />
     76        <field field_ref="BIN23" level="2" />
     77        <field field_ref="BIN24" level="2" />
     78        <field field_ref="BIN25" level="2" />
     79        <field field_ref="BIN26" level="2" />
     80        <field field_ref="BIN27" level="2" />
     81        <field field_ref="BIN28" level="2" />
     82        <field field_ref="BIN29" level="2" />
     83        <field field_ref="BIN30" level="2" />
     84        <field field_ref="BIN31" level="2" />
     85        <field field_ref="BIN32" level="2" />
     86        <field field_ref="BIN33" level="2" />
     87        <field field_ref="BIN34" level="2" />
     88        <field field_ref="BIN35" level="2" />
     89        <field field_ref="BIN36" level="2" />
     90      </field_group>
     91     
     92    </file>
     93  </file_group>
     94</file_definition>
    495
    5             <field_group group_ref="fields_strataer_3D" operation="average" freq_op="1ts" level="1" />
    6             <field_group group_ref="fields_strataer_2D" operation="average" freq_op="1ts" level="1" />
    7             <field_group group_ref="fields_strataer_trac_3D" operation="average" freq_op="1ts" level="1" />
    8             <field_group group_ref="fields_strataer_trac_2D" operation="average" freq_op="1ts" level="1" />
    9            
    10         </file>
    11     </file_group>
    12 </file_definition>

  • LMDZ6/trunk/DefLists/file_def_histmth_lmdz.xml

    r4298 r4601  
    675675                <field field_ref="cool_volc" level="1" />
    676676                <field field_ref="heat_volc" level="1" />
     677                <field field_ref="dqch4"     level="10" />
    677678            </field_group>
    678679

  • LMDZ6/trunk/DefLists/file_def_histstrataer_lmdz.xml

    r3507 r4601  
    11<!-- $Id$ -->
    22<file_definition>
    3     <file_group id="defile">
    4         <file id="histstrataer" name="histstrataer" output_freq="1mo" output_level="_AUTO_" enabled="_AUTO_" compression_level="2">
    5 
    6             <field_group group_ref="remap_1d">
    7               <field_group group_ref="fields_strataer_3D" grid_ref="grid_out_presnivs" level="1" />
    8               <field_group group_ref="fields_strataer_2D" grid_ref="grid_out"  level="1" />
    9               <field_group group_ref="fields_strataer_trac_3D"  grid_ref="grid_out_presnivs" level="1" />
    10               <field_group group_ref="fields_strataer_trac_2D" grid_ref="grid_out" level="1" />
    11            </field_group>
    12            
    13         </file>
    14     </file_group>
     3  <file_group id="defile">
     4    <file id="histstrataer" name="histstrataer" output_freq="1mo" output_level="_AUTO_" enabled="_AUTO_" compression_level="2">
     5     
     6      <field_group group_ref="remap_1mo" operation="average">
     7       
     8        <field_group grid_ref="grid_out"  level="3">
     9          <field field_ref="OD550_strat_only" level="1" />
     10          <field field_ref="OD1020_strat_only" level="1" />
     11          <field field_ref="surf_PM25_sulf" level="1" />
     12          <field field_ref="budg_dep_dry_ocs" level="3" />
     13          <field field_ref="budg_dep_wet_ocs" level="3" />
     14          <field field_ref="budg_dep_dry_so2" level="3" />
     15          <field field_ref="budg_dep_wet_so2" level="3" />
     16          <field field_ref="budg_dep_dry_h2so4" level="1" />
     17          <field field_ref="budg_dep_wet_h2so4" level="1" />
     18          <field field_ref="budg_dep_dry_part" level="1" />
     19          <field field_ref="budg_dep_wet_part" level="1" />
     20          <field field_ref="budg_emi_ocs" level="1" />
     21          <field field_ref="budg_emi_so2" level="1" />
     22          <field field_ref="budg_emi_h2so4" level="1" />
     23          <field field_ref="budg_emi_part" level="1" />
     24          <field field_ref="budg_ocs_to_so2" level="1" />
     25          <field field_ref="budg_so2_to_h2so4" level="1" />
     26          <field field_ref="budg_h2so4_to_part" level="1" />
     27          <field field_ref="budg_sed_part" level="1" />
     28          <field field_ref="p_tropopause" level="1" />
     29          <field field_ref="z_tropopause" level="1" />
     30          <field field_ref="t_tropopause" level="3" />
     31        </field_group>
     32       
     33        <field_group grid_ref="grid_out_presnivs" level="5">
     34          <field field_ref="ext_strat_550"  level="1" />
     35          <field field_ref="ext_strat_1020" level="1" />
     36          <field field_ref="budg_3D_nucl" level="1" />
     37          <field field_ref="budg_3D_cond_evap" level="1" />
     38          <field field_ref="budg_3D_ocs_to_so2" level="1" />
     39          <field field_ref="budg_3D_so2_to_h2so4" level="1" />
     40          <field field_ref="budg_3D_backgr_ocs" level="1" />
     41          <field field_ref="budg_3D_backgr_so2" level="1" />
     42          <field field_ref="R2SO4" level="1" />
     43          <field field_ref="OCS_lifetime" level="1" />
     44          <field field_ref="SO2_lifetime" level="1" />
     45          <field field_ref="vsed_aer" level="2" />
     46          <field field_ref="f_r_wet" level="1" />
     47          <field field_ref="mass" level="2" />
     48          <field field_ref="temp" level="2" />
     49          <field field_ref="pres" level="2" />
     50          <field field_ref="h2o"  level="1" />
     51        </field_group>
     52       
     53        <field_group grid_ref="grid_out_presnivs" level="5">
     54          <field field_ref="BIN01" level="1" />
     55          <field field_ref="BIN02" level="1" />
     56          <field field_ref="BIN03" level="1" />
     57          <field field_ref="BIN04" level="1" />
     58          <field field_ref="BIN05" level="1" />
     59          <field field_ref="BIN06" level="1" />
     60          <field field_ref="BIN07" level="1" />
     61          <field field_ref="BIN08" level="1" />
     62          <field field_ref="BIN09" level="1" />
     63          <field field_ref="BIN10" level="1" />
     64          <field field_ref="BIN11" level="1" />
     65          <field field_ref="BIN12" level="1" />
     66          <field field_ref="BIN13" level="1" />
     67          <field field_ref="BIN14" level="1" />
     68          <field field_ref="BIN15" level="1" />
     69          <field field_ref="BIN16" level="1" />
     70          <field field_ref="BIN17" level="1" />
     71          <field field_ref="BIN18" level="1" />
     72          <field field_ref="BIN19" level="1" />
     73          <field field_ref="BIN20" level="1" />
     74          <field field_ref="BIN21" level="1" />
     75          <field field_ref="BIN22" level="1" />
     76          <field field_ref="BIN23" level="1" />
     77          <field field_ref="BIN24" level="1" />
     78          <field field_ref="BIN25" level="1" />
     79          <field field_ref="BIN26" level="1" />
     80          <field field_ref="BIN27" level="1" />
     81          <field field_ref="BIN28" level="1" />
     82          <field field_ref="BIN29" level="1" />
     83          <field field_ref="BIN30" level="1" />
     84          <field field_ref="BIN31" level="1" />
     85          <field field_ref="BIN32" level="1" />
     86          <field field_ref="BIN33" level="1" />
     87          <field field_ref="BIN34" level="1" />
     88          <field field_ref="BIN35" level="1" />
     89          <field field_ref="BIN36" level="1" />
     90          <field field_ref="GASOCS"  level="1" />
     91          <field field_ref="GASSO2"  level="1" />
     92          <field field_ref="GASH2SO4"  level="1" />
     93        </field_group>
     94       
     95      </field_group>
     96    </file>
     97  </file_group>
    1598</file_definition>

  • LMDZ6/trunk/libf/phylmd/StratAer/aerophys.F90

    r3526 r4601  
    1717  REAL, PARAMETER                        :: mSatom=32.06*1.66E-27    ! Mass of a S atom [kg]
    1818  REAL, PARAMETER                        :: mOCSmol=60.07*1.66E-27   ! Mass of an OCS molecule [kg]
     19  REAL, PARAMETER                        :: mClatom=35.45*1.66E-27   ! Mass of an Cl atom [kg]
     20  REAL, PARAMETER                        :: mHClmol=36.46*1.66E-27   ! Mass of an HCl molecule [kg]
     21  REAL, PARAMETER                        :: mBratom=79.90*1.66E-27   ! Mass of an Br atom [kg]
     22  REAL, PARAMETER                        :: mHBrmol=80.92*1.66E-27   ! Mass of an HBr molecule [kg]
     23  REAL, PARAMETER                        :: mNOmol=30.01*1.66E-27    ! Mass of an NO molecule [kg]
     24  REAL, PARAMETER                        :: mNO2mol=46.01*1.66E-27   ! Mass of an NO2 molecule [kg]
     25  REAL, PARAMETER                        :: mNatome=14.0067*1.66E-27 ! Mass of an N atome [kg]
    1926!
    2027END MODULE aerophys

  • LMDZ6/trunk/libf/phylmd/StratAer/interp_sulf_input.F90

    r3677 r4601  
    1212  USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
    1313  USE phys_local_var_mod, ONLY : budg_3D_backgr_ocs, budg_3D_backgr_so2
    14   USE phys_local_var_mod, ONLY : OCS_lifetime, SO2_lifetime
     14  USE phys_local_var_mod, ONLY : OCS_lifetime, SO2_lifetime, H2SO4_lifetime, O3_clim
    1515  USE mod_phys_lmdz_para
    1616  USE dimphy
    … …  
    1919  USE aerophys
    2020  USE YOMCST
     21  USE strataer_local_var_mod, ONLY : flag_newclim_file,flag_verbose_strataer
    2122
    2223  IMPLICIT NONE
    … …  
    6263  REAL OCS_clim_glo(klon_glo,klev)
    6364  REAL SO2_clim_glo(klon_glo,klev)
     65 
     66  REAL, ALLOCATABLE :: O3_clim_in(:, :, :, :)
     67  REAL, ALLOCATABLE :: O3_clim_mth(:, :, :)
     68  REAL, ALLOCATABLE :: O3_clim_tmp(:, :)
     69  REAL O3_clim_glo(klon_glo,klev)
     70 
     71  ! fixed climatos
    6472  REAL, ALLOCATABLE :: OCS_lifetime_in(:, :, :, :)
    6573  REAL, ALLOCATABLE :: SO2_lifetime_in(:, :, :, :)
    … …  
    7078  REAL OCS_lifetime_glo(klon_glo,klev)
    7179  REAL SO2_lifetime_glo(klon_glo,klev)
     80 
     81  REAL, ALLOCATABLE :: H2SO4_lifetime_in(:, :, :, :)
     82  REAL, ALLOCATABLE :: H2SO4_lifetime_mth(:, :, :)
     83  REAL, ALLOCATABLE :: H2SO4_lifetime_tmp(:, :)
     84
     85  REAL H2SO4_lifetime_glo(klon_glo,klev)
    7286!
    7387  REAL, ALLOCATABLE, SAVE :: OCS_clim(:,:)
    … …  
    7993
    8094! For NetCDF:
    81   INTEGER ncid_in  ! IDs for input files
    82   INTEGER varid, ncerr
     95  INTEGER             :: ncid_in  ! IDs for input files
     96  INTEGER             :: varid, ncerr
     97  CHARACTER (LEN=3)   :: nc_lat, nc_lon
     98  CHARACTER (LEN=256) :: nc_fname
    8399   
    84100  INTEGER, PARAMETER :: lev_input=17
    … …  
    103119  IF (is_mpi_root.AND.is_omp_root) THEN
    104120
    105 !--reading emission files
    106     CALL nf95_open("ocs_so2_annual_lmdz.nc", nf90_nowrite, ncid_in)
    107 
     121    !--init ncdf variables
     122    IF(flag_newclim_file) THEN
     123      nc_fname = "ocs_so2_h2so4_annual_lmdz.nc"
     124      nc_lat   = "LAT"
     125      nc_lon   = "LON"
     126    ELSE
     127      ! old file for retro compatibility
     128      nc_fname = "ocs_so2_annual_lmdz.nc"
     129      nc_lat   = "lat"
     130      nc_lon   = "lon"
     131    ENDIF
     132
     133    !--reading emission files
     134    CALL nf95_open(nc_fname, nf90_nowrite, ncid_in)
     135    IF(flag_verbose_strataer) print *, "Open file ", nc_fname
     136   
    108137    CALL nf95_inq_varid(ncid_in, "LEV", varid)
    109138    CALL nf95_gw_var(ncid_in, varid, lev)
    110139    n_lev = size(lev)
    111 
    112     CALL nf95_inq_varid(ncid_in, "lat", varid)
     140   
     141    CALL nf95_inq_varid(ncid_in, nc_lat, varid)
    113142    CALL nf95_gw_var(ncid_in, varid, latitude)
    114143    n_lat = size(latitude)
    115 
    116     CALL nf95_inq_varid(ncid_in, "lon", varid)
     144   
     145    CALL nf95_inq_varid(ncid_in, nc_lon, varid)
    117146    CALL nf95_gw_var(ncid_in, varid, longitude)
    118147    n_lon = size(longitude)
    119 
     148   
    120149    CALL nf95_inq_varid(ncid_in, "TIME", varid)
    121150    CALL nf95_gw_var(ncid_in, varid, time)
    122151    n_mth = size(time)
    123 
     152   
    124153    IF (.NOT.ALLOCATED(OCS_clim_in))     ALLOCATE(OCS_clim_in(n_lon, n_lat, n_lev, n_mth))
    125154    IF (.NOT.ALLOCATED(SO2_clim_in))     ALLOCATE(SO2_clim_in(n_lon, n_lat, n_lev, n_mth))
     155    IF (.NOT.ALLOCATED(O3_clim_in))      ALLOCATE(O3_clim_in(n_lon, n_lat, n_lev, n_mth))
     156   
    126157    IF (.NOT.ALLOCATED(OCS_lifetime_in)) ALLOCATE(OCS_lifetime_in(n_lon, n_lat, n_lev, n_mth))
    127158    IF (.NOT.ALLOCATED(SO2_lifetime_in)) ALLOCATE(SO2_lifetime_in(n_lon, n_lat, n_lev, n_mth))
    128 
     159    IF (.NOT.ALLOCATED(H2SO4_lifetime_in)) ALLOCATE(H2SO4_lifetime_in(n_lon, n_lat, n_lev, n_mth))
     160   
    129161    CALL nf95_inq_varid(ncid_in, "OCS", varid)
    130162    ncerr = nf90_get_var(ncid_in, varid, OCS_clim_in)
    131     print *,'code erreur OCS=', ncerr, varid
     163    IF(flag_verbose_strataer) print *,'code erreur OCS=', ncerr, varid
    132164
    133165    CALL nf95_inq_varid(ncid_in, "SO2", varid)
    134166    ncerr = nf90_get_var(ncid_in, varid, SO2_clim_in)
    135     print *,'code erreur SO2=', ncerr, varid
     167    IF(flag_verbose_strataer) print *,'code erreur SO2=', ncerr, varid
    136168
    137169    CALL nf95_inq_varid(ncid_in, "OCS_LIFET", varid)
    138170    ncerr = nf90_get_var(ncid_in, varid, OCS_lifetime_in)
    139     print *,'code erreur OCS_lifetime_in=', ncerr, varid
     171    IF(flag_verbose_strataer) print *,'code erreur OCS_lifetime_in=', ncerr, varid
    140172
    141173    CALL nf95_inq_varid(ncid_in, "SO2_LIFET", varid)
    142174    ncerr = nf90_get_var(ncid_in, varid, SO2_lifetime_in)
    143     print *,'code erreur SO2_lifetime_in=', ncerr, varid
    144 
     175    IF(flag_verbose_strataer) print *,'code erreur SO2_lifetime_in=', ncerr, varid
     176   
     177    IF(flag_newclim_file) THEN
     178       CALL nf95_inq_varid(ncid_in, "O3", varid)
     179       ncerr = nf90_get_var(ncid_in, varid, O3_clim_in)
     180       IF(flag_verbose_strataer) print *,'code erreur O3=', ncerr, varid
     181       
     182       CALL nf95_inq_varid(ncid_in, "H2SO4_LIFET", varid)
     183       ncerr = nf90_get_var(ncid_in, varid, H2SO4_lifetime_in)
     184       IF(flag_verbose_strataer) print *,'code erreur H2SO4_lifetime_in=', ncerr, varid
     185    ENDIF
     186   
    145187    CALL nf95_close(ncid_in)
    146188
    … …  
    149191    IF (.NOT.ALLOCATED(OCS_clim_tmp)) ALLOCATE(OCS_clim_tmp(klon_glo, n_lev))
    150192    IF (.NOT.ALLOCATED(SO2_clim_tmp)) ALLOCATE(SO2_clim_tmp(klon_glo, n_lev))
     193    IF (.NOT.ALLOCATED(O3_clim_mth))  ALLOCATE(O3_clim_mth(n_lon, n_lat, n_lev))
     194    IF (.NOT.ALLOCATED(O3_clim_tmp))  ALLOCATE(O3_clim_tmp(klon_glo, n_lev))
     195   
    151196    IF (.NOT.ALLOCATED(OCS_lifetime_mth)) ALLOCATE(OCS_lifetime_mth(n_lon, n_lat, n_lev))
    152197    IF (.NOT.ALLOCATED(SO2_lifetime_mth)) ALLOCATE(SO2_lifetime_mth(n_lon, n_lat, n_lev))
    153198    IF (.NOT.ALLOCATED(OCS_lifetime_tmp)) ALLOCATE(OCS_lifetime_tmp(klon_glo, n_lev))
    154199    IF (.NOT.ALLOCATED(SO2_lifetime_tmp)) ALLOCATE(SO2_lifetime_tmp(klon_glo, n_lev))
    155 
     200    IF (.NOT.ALLOCATED(H2SO4_lifetime_mth)) ALLOCATE(H2SO4_lifetime_mth(n_lon, n_lat, n_lev))
     201    IF (.NOT.ALLOCATED(H2SO4_lifetime_tmp)) ALLOCATE(H2SO4_lifetime_tmp(klon_glo, n_lev))
     202   
    156203!---select the correct month, undo multiplication with 1.e12 (precision reasons)
    157204!---correct latitudinal order and convert input from volume mixing ratio to mass mixing ratio
    … …  
    161208      SO2_lifetime_mth(:,j,:) = SO2_lifetime_in(:,n_lat+1-j,:,mth_cur)
    162209      OCS_lifetime_mth(:,j,:) = OCS_lifetime_in(:,n_lat+1-j,:,mth_cur)
     210     
     211      ! O3 from 2d model is not tracer, in VMR
     212      IF(flag_newclim_file) THEN
     213         H2SO4_lifetime_mth(:,j,:) = H2SO4_lifetime_in(:,n_lat+1-j,:,mth_cur)
     214         ! new input files
     215         O3_clim_mth(:,j,:) = 1.e-6*O3_clim_in(:,n_lat+1-j,:,mth_cur)
     216      ELSE
     217         H2SO4_lifetime_mth(:,j,:) = 1.e-6
     218         O3_clim_mth(:,j,:) = 1.e-6
     219      ENDIF
    163220    ENDDO
    164221
    … …  
    166223    CALL grid2dTo1d_glo(OCS_clim_mth,OCS_clim_tmp)
    167224    CALL grid2dTo1d_glo(SO2_clim_mth,SO2_clim_tmp)
     225    CALL grid2dTo1d_glo(O3_clim_mth,O3_clim_tmp)
     226   
    168227    CALL grid2dTo1d_glo(OCS_lifetime_mth,OCS_lifetime_tmp)
    169228    CALL grid2dTo1d_glo(SO2_lifetime_mth,SO2_lifetime_tmp)
     229    CALL grid2dTo1d_glo(H2SO4_lifetime_mth,H2SO4_lifetime_tmp)
    170230
    171231!--set lifetime to very high value in uninsolated areas
    172232    DO i=1, klon_glo
    173233      DO kk=1, n_lev
     234!--added tests on VMR of species for realism
     235        IF (OCS_clim_tmp(i,kk).LT.1.e-20) THEN
     236          OCS_clim_tmp(i,kk)=1.0e-20
     237        ENDIF
     238        IF (SO2_clim_tmp(i,kk).LT.1.e-20) THEN
     239          SO2_clim_tmp(i,kk)=1.0e-20
     240        ENDIF
     241        !       50 ppbv min for O3
     242        IF (O3_clim_tmp(i,kk).LT.50.0e-9) THEN
     243           O3_clim_tmp(i,kk)=50.0e-9
     244        ENDIF
     245       
    174246        IF (OCS_lifetime_tmp(i,kk)==0.0) THEN
    175247          OCS_lifetime_tmp(i,kk)=1.0e12
    … …  
    178250          SO2_lifetime_tmp(i,kk)=1.0e12
    179251        ENDIF
     252       IF (H2SO4_lifetime_tmp(i,kk)==0.0) THEN
     253          H2SO4_lifetime_tmp(i,kk)=1.0e12
     254       ENDIF
    180255      ENDDO
    181256    ENDDO
    … …  
    186261     OCS_lifetime_glo(i,k)=0.0
    187262     SO2_lifetime_glo(i,k)=0.0
     263     O3_clim_glo(i,k)=0.0
     264     
    188265     OCS_clim_glo(i,k)=0.0
    189266     SO2_clim_glo(i,k)=0.0
     267     H2SO4_lifetime_glo(i,k)=0.0
     268     
    190269     DO kk=1, n_lev
    191270      OCS_lifetime_glo(i,k)=OCS_lifetime_glo(i,k)+ &
    … …  
    195274           MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
    196275           *SO2_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
     276      IF(flag_newclim_file) THEN
     277         H2SO4_lifetime_glo(i,k)=H2SO4_lifetime_glo(i,k)+ &
     278              MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) &
     279              -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
     280              *H2SO4_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
     281      ENDIF
     282     
    197283      OCS_clim_glo(i,k)=OCS_clim_glo(i,k)+ &
    198284           MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
    … …  
    201287           MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
    202288           *SO2_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
     289      IF(flag_newclim_file) THEN
     290         O3_clim_glo(i,k)=O3_clim_glo(i,k)+ &
     291              MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) &
     292              -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) &
     293              *O3_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1))
     294      ENDIF
    203295      ENDDO
    204296    ENDDO
    … …  
    213305  CALL scatter(OCS_clim_glo, OCS_clim)
    214306  CALL scatter(SO2_clim_glo, SO2_clim)
     307  CALL scatter(O3_clim_glo, O3_clim)
    215308  CALL scatter(OCS_lifetime_glo, OCS_lifetime)
    216309  CALL scatter(SO2_lifetime_glo, SO2_lifetime)
    217 
     310  CALL scatter(H2SO4_lifetime_glo, H2SO4_lifetime)
     311 
    218312  IF (is_mpi_root.AND.is_omp_root) THEN
    219313!
    220     DEALLOCATE(OCS_clim_in,SO2_clim_in)
    221     DEALLOCATE(OCS_clim_mth,SO2_clim_mth)
    222     DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp)
    223     DEALLOCATE(OCS_lifetime_in,SO2_lifetime_in)
    224     DEALLOCATE(OCS_lifetime_mth,SO2_lifetime_mth)
    225     DEALLOCATE(OCS_lifetime_tmp,SO2_lifetime_tmp)
     314    DEALLOCATE(OCS_clim_in,SO2_clim_in,O3_clim_in)
     315    DEALLOCATE(OCS_clim_mth,SO2_clim_mth,O3_clim_mth)
     316    DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp,O3_clim_tmp)
     317    DEALLOCATE(OCS_lifetime_in, SO2_lifetime_in, H2SO4_lifetime_in)
     318    DEALLOCATE(OCS_lifetime_mth, SO2_lifetime_mth, H2SO4_lifetime_mth)
     319    DEALLOCATE(OCS_lifetime_tmp, SO2_lifetime_tmp, H2SO4_lifetime_tmp)
    226320!
    227321  ENDIF !--is_mpi_root and is_omp_root

  • LMDZ6/trunk/libf/phylmd/StratAer/micphy_tstep.F90

    r4293 r4601  
    1414  USE YOMCST, ONLY : RPI, RD, RG
    1515  USE print_control_mod, ONLY: lunout
    16   USE strataer_mod
     16  USE strataer_local_var_mod
    1717 
    1818  IMPLICIT NONE

  • LMDZ6/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90

    r4293 r4601  
    1111  USE infotrac_phy
    1212  USE YOMCST, ONLY : RPI, RD, RMD, RKBOL, RNAVO
    13 
     13  USE strataer_local_var_mod, ONLY : flag_new_nucl
     14 
    1415  IMPLICIT NONE
    1516
    1617  ! input variables
    17   LOGICAL, PARAMETER :: flag_new_nucl=.TRUE.
    1818  REAL rhoa    ! H2SO4 number density [molecules/cm3]
    1919  REAL t_seri  ! temperature (K)
    … …  
    170170
    171171  IF (t < 190.15) THEN
    172 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature < 190.15 K, using 190.15 K (T=',t,'K)'
    173172     t=190.15
    174173  ENDIF
    175174
    176175  IF (t > 300.15) THEN
    177 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature > 300.15 K, using 300.15 K'
    178176     t=300.15
    179177  ENDIF
    180178
    181179  IF (rh < 0.0001) THEN
    182 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water < 0.0001, using 0.0001'
    183180     rh=0.0001
    184181  ENDIF
    185182
    186183  IF (rh > 1.0) THEN
    187 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water > 1 using 1'
    188 !     print *, 'rh=',rh
    189184     rh=1.0
    190185  ENDIF
    191186
    192187  IF (rhoa < 1.e4) THEN
    193 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration < 1e4 1/cm3, using 1e4 1/cm3'
    194188     rhoa=1.e4
    195189  ENDIF
    196190
    197191  IF (rhoa > 1.e11) THEN
    198 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration > 1e11 1/cm3, using 1e11 1/cm3'
    199192     rhoa=1.e11
    200193  ENDIF
    … …  
    281274
    282275  IF (jnuc < 1.E-7) THEN
    283 !     print *,'Warning (ilon=',ilon,'ilev=',ilev'): nucleation rate < 1E-7/cm3s, using 0.0/cm3s,'
    284276     jnuc=0.0
    285277  ENDIF
    286278
    287279  IF (jnuc > 1.e10) THEN
    288 !     print *,'Warning (ilon=',ilon,'ilev=',ilev, &
    289 !        & '): nucleation rate > 1e10/cm3s, using 1e10/cm3s'
    290280     jnuc=1.e10
    291281  ENDIF
    … …  
    425415  !Temperature bounds
    426416  IF (t.LE.165.) THEN
    427 !     print *,'Warning: temperature < 165.0 K, using 165.0 K in neutral nucleation calculation'
    428417     tln=165.0
    429418  ENDIF
    430419  IF (t.LE.195.) THEN
    431 !     print *,'Warning: temperature < 195.0 K, using 195.0 K in ion-induced nucleation calculation'
    432420     tli=195.0
    433421  ENDIF
    434422  IF (t.GE.400.) THEN
    435 !     print *,'Warning: temperature > 400. K, using 400. K in nucleation calculations'
    436423     tln=400.
    437424     tli=400.
    … …  
    440427  ! Saturation ratio bounds
    441428  IF (satrat.LT.1.E-7) THEN
    442 !     print *,'Warning: saturation ratio of water < 1.E-7, using 1.E-7 in ion-induced nucleation calculation'
    443429     satratli=1.E-7
    444430  ENDIF
    445431  IF (satrat.LT.1.E-5) THEN
    446 !     print *,'Warning: saturation ratio of water < 1.E-5, using 1.E-5 in neutral nucleation calculation'
    447432     satratln=1.E-5
    448433  ENDIF
    449434  IF (satrat.GT.0.95) THEN
    450 !     print *,'Warning: saturation ratio of water > 0.95, using 0.95 in ion-induced nucleation calculation'
    451435     satratli=0.95
    452436  ENDIF
    453437  IF (satrat.GT.1.0) THEN
    454 !     print *,'Warning: saturation ratio of water > 1 using 1 in neutral nucleation calculation'
    455438     satratln=1.0
    456439  ENDIF
    … …  
    458441  ! Sulfuric acid concentration bounds
    459442  IF (rhoa.LE.1.E4) THEN
    460 !     print *,'Warning: sulfuric acid < 1e4 1/cm3, using 1e4 1/cm3 in nucleation calculation'
    461443     rhoaln=1.E4
    462444     rhoali=1.E4
    463445  ENDIF
    464446  IF (rhoa.GT.1.E13) THEN
    465 !     print *,'Warning: sulfuric acid > 1e13 1/cm3, using 1e13 1/cm3 in neutral nucleation calculation'
    466447     rhoaln=1.E13
    467448  ENDIF
    468449  IF (rhoa.GT.1.E16) THEN
    469 !     print *,'Warning: sulfuric acid concentration > 1e16 1/cm3, using 1e16 1/cm3 in ion-induced nucleation calculation'
    470450     rhoali=1.E16
    471451  ENDIF
    … …  
    521501     ! Dimer formation rate
    522502     jnuc_n=1.E6*(2.*0.3E-9)**2.*SQRT(8.*RPI*RKBOL*(1./mH2SO4mol+1./mH2SO4mol))/2.*SQRT(t)*rhoa**2.
    523 !     jnuc_n=1.E6*(2.*0.3E-9)**2.*SQRT(8.*3.141593*1.38E-23*(1./(1.661E-27*98.07)+1./(1.661E-27*98.07)))/2.*SQRT(t)*rhoa**2.
    524503     ntot_n=1. !set to 1
    525504     na_n=1.   ! The critical cluster contains one molecules but the produced cluster contains 2 molecules
    … …  
    604583     na_n=x_n*ntot_n
    605584     IF (na_n .LT. 1.) THEN
    606 !        print *, 'Warning: number of acid molecules < 1 in nucleation regime, setting na_n=1'
    607585        na_n=1.0
    608586     ENDIF
    … …  
    664642     
    665643     IF (kinetic_i) THEN   
    666 !        jnuc_i1=1.0E6*(0.3E-9 + 0.487E-9)**2.*SQRT(8.*3.141593*1.38E-23*(1./(1.661E-27*98.07)+1./(1.661E-27*98.07)))*  &
    667644        jnuc_i1=1.0E6*(0.3E-9 + 0.487E-9)**2.*SQRT(8.*RPI*RKBOL*(1./mH2SO4mol+1./mH2SO4mol))*  &
    668645             &  SQRT(tli)*rhoali !1/cm3s 
    … …  
    816793        na_i=x_i*ntot_i
    817794        IF (na_i .LT. 1.) THEN
    818 !           print *, 'Warning: number of acid molecules < 1 in nucleation regime, setting na_n=1'
    819795           na_n=1.0
    820796        ENDIF

  • LMDZ6/trunk/libf/phylmd/StratAer/ocs_to_so2.F90

    r3677 r4601  
    88  USE infotrac_phy
    99  USE YOMCST, ONLY : RG
    10   USE phys_local_var_mod, ONLY : OCS_lifetime, budg_3D_ocs_to_so2, budg_ocs_to_so2
     10  USE phys_local_var_mod, ONLY : OCS_lifetime, budg_3D_ocs_to_so2, budg_ocs_to_so2
     11  USE strataer_local_var_mod, ONLY : flag_min_rreduce
    1112
    1213  IMPLICIT NONE
    … …  
    2324  ! local variables
    2425  INTEGER                                       :: i,j,k,nb,ilon,ilev
    25 
     26  REAL                                          :: rreduce
     27 
    2628!--convert OCS to SO2
    2729  budg_3D_ocs_to_so2(:,:)=0.0
    … …  
    2931
    3032  DO ilon=1, klon
    31   DO ilev=1, klev
    32   !only in the stratosphere
    33   IF (is_strato(ilon,ilev)) THEN
    34     IF (OCS_lifetime(ilon,ilev).GT.0.0) THEN
    35       budg_3D_ocs_to_so2(ilon,ilev)=tr_seri(ilon,ilev,id_OCS_strat)*(1.0-exp(-pdtphys/OCS_lifetime(ilon,ilev)))
    36     ENDIF
    37     tr_seri(ilon,ilev,id_OCS_strat)=tr_seri(ilon,ilev,id_OCS_strat) - budg_3D_ocs_to_so2(ilon,ilev)
    38     tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + mSO2mol/mOCSmol*budg_3D_ocs_to_so2(ilon,ilev)
    39     !convert budget from kg(OCS)/kgA to kg(S)/m2/layer/s for saving as diagnostic
    40     budg_3D_ocs_to_so2(ilon,ilev)=budg_3D_ocs_to_so2(ilon,ilev)*mSatom/mOCSmol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
    41     budg_ocs_to_so2(ilon)=budg_ocs_to_so2(ilon)+budg_3D_ocs_to_so2(ilon,ilev)
    42   ENDIF
    43   ENDDO
     33     DO ilev=1, klev
     34        !only in the stratosphere
     35        IF (is_strato(ilon,ilev)) THEN
     36          IF (OCS_lifetime(ilon,ilev).GT.0.0.AND.OCS_lifetime(ilon,ilev).LT.1.E10) THEN
     37            rreduce = OCS_lifetime(ilon,ilev)
     38            ! Check lifetime rreduce < timestep*3 (such as H2SO4 loss > 0.28*H2SO4) with exp(-1/3)=0.72
     39            IF(flag_min_rreduce) THEN
     40               IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys
     41            ENDIF
     42            budg_3D_ocs_to_so2(ilon,ilev)=tr_seri(ilon,ilev,id_OCS_strat)*(1.0-exp(-pdtphys/rreduce))
     43            tr_seri(ilon,ilev,id_OCS_strat)=tr_seri(ilon,ilev,id_OCS_strat) - budg_3D_ocs_to_so2(ilon,ilev)
     44            tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + mSO2mol/mOCSmol*budg_3D_ocs_to_so2(ilon,ilev)
     45            !convert budget from kg(OCS)/kgA to kg(S)/m2/layer/s for saving as diagnostic
     46            budg_3D_ocs_to_so2(ilon,ilev)=budg_3D_ocs_to_so2(ilon,ilev)*mSatom/mOCSmol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
     47            budg_ocs_to_so2(ilon)=budg_ocs_to_so2(ilon)+budg_3D_ocs_to_so2(ilon,ilev)
     48         ENDIF
     49         ! END IF(OCS_lifetime(ilon,ilev).GT.0.0.AND.OCS_lifetime(ilon,ilev).LT.1.E10)
     50      ENDIF
     51      ! END IF(is_strato(ilon,ilev))
     52   ENDDO
    4453  ENDDO
    4554

  • LMDZ6/trunk/libf/phylmd/StratAer/so2_to_h2so4.F90

    r3677 r4601  
    77  USE aerophys
    88  USE infotrac_phy
    9   USE YOMCST, ONLY : RG
    10   USE phys_local_var_mod, ONLY : SO2_lifetime, budg_3D_so2_to_h2so4, budg_so2_to_h2so4
    11 
     9  USE YOMCST, ONLY : RG, RD
     10  ! lifetime (sec) et O3_clim (VMR)
     11  USE phys_local_var_mod, ONLY : SO2_lifetime, H2SO4_lifetime, O3_clim, budg_3D_so2_to_h2so4, budg_so2_to_h2so4
     12  USE strataer_local_var_mod, ONLY : flag_OH_reduced, flag_H2SO4_photolysis, flag_min_rreduce
     13 
    1214  IMPLICIT NONE
    1315
    … …  
    2426  ! local variables in coagulation routine
    2527  INTEGER                                       :: i,j,k,nb,ilon,ilev
    26 
     28  REAL                                          :: dummyso4toso2,rkho2o3,rkoho3,rkso2oh
     29  REAL                                          :: rmairdens,rrak0,rrak1,rrate,rreduce
     30 
    2731!--convert SO2 to H2SO4
    2832  budg_3D_so2_to_h2so4(:,:)=0.0
    … …  
    3034
    3135  DO ilon=1, klon
    32   DO ilev=1, klev
    33   !only in the stratosphere
    34   IF (is_strato(ilon,ilev)) THEN
    35     IF (SO2_lifetime(ilon,ilev).GT.0.0) THEN
    36       budg_3D_so2_to_h2so4(ilon,ilev)=tr_seri(ilon,ilev,id_SO2_strat)*(1.0-exp(-pdtphys/SO2_lifetime(ilon,ilev)))
    37     ENDIF
    38     tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - budg_3D_so2_to_h2so4(ilon,ilev)
    39     tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) + mH2SO4mol/mSO2mol*budg_3D_so2_to_h2so4(ilon,ilev)
    40     !convert budget from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic
    41     budg_3D_so2_to_h2so4(ilon,ilev)=budg_3D_so2_to_h2so4(ilon,ilev)*mSatom/mSO2mol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
    42     budg_so2_to_h2so4(ilon)=budg_so2_to_h2so4(ilon)+budg_3D_so2_to_h2so4(ilon,ilev)
    43   ENDIF
    44   ENDDO
    45   ENDDO
    46 
     36     DO ilev=1, klev
     37     !only in the stratosphere
     38        IF (is_strato(ilon,ilev)) THEN
     39           IF (SO2_lifetime(ilon,ilev).GT.0.0 .AND. SO2_lifetime(ilon,ilev).LT.1.E10) THEN
     40              IF (flag_OH_reduced) THEN
     41                 !--convert SO2 to H2SO4 (slimane)
     42                 ! Reduce OH because of SO2
     43                 ! d[OH]/dt = k1[HO2][O] +k2[HO2][O3] +k3[HO2][NO]
     44                 ! − (k4[O] +k5[O3] +k6[CO] +ks[SO2])[OH]
     45                 ! In steady-state: d[OH]/dt = 0, so the ratio [OH]/[HO2] =
     46                 ! (k1[O] +k2[O3] +k3[NO]) / (k4[O] +k5[O3] +k6[CO] +ks[SO2])
     47                 ! In the lower and middle stratosphere,
     48                 ! [OH]/[HO2] ~(k2[O3] +k3[NO]) / (k5[O3] +k6[CO] +ks[SO2])
     49                 ! CO is only important near the tropopause. NO should not dominate over O3.
     50                 ! So when consideringthe big terms: [OH]/[HO2] ~ ~ k2[O3] / (k5[O3] +ks[SO2])
     51                 ! For c: control run c (no effect of SO2), [OHc]/[HO2c] ~ ~ k2[O3] / k5[O3] = 1/ Rc
     52                 ! For p: perturbed run (high SO2), [OHp]/[HO2p] ~ ~ k2[O3] / (k5[O3] +ks[SO2]) =1/Rp
     53                 ! If HOx = OH + HO2 = constante,
     54                 ! OHc + HO2c = OHp + HO2p
     55                 ! OHc (1+Rc) = OHp (1+Rp)
     56                 ! OHp = OHc (1+Rc) / (1+Rp)
     57                 ! 1+Rc = (k2[O3] +k5[O3] ) / k2[O3]
     58                 ! 1+Rp = (k2[O3] +k5[O3] +ks[SO2]) / k2[O3]
     59                 ! (1+Rc) / (1+Rp)= (k2[O3] +k5[O3] )/ (k2[O3] +k5[O3] +ks[SO2])
     60                 ! OHp = OHc * (k(HO2+OH)*[O3] +k(OH+O3)[O3] ) /
     61                 ! (k(HO2+OH)*[O3] +k(OH+O3)[O3] +k(SO2+OH)[SO2])
     62                 ! Final: OHp = OHc * (k(HO2+O3)*[O3] +k(OH+O3)[O3] ) /
     63                 !  (k(HO2+O3)*[O3] +k(OH+O3)[O3] +k(SO2+OH)[SO2])
     64                 
     65                 !   latest jpl 2020
     66                 rkho2o3 = 1.0e-14*exp( -490./t_seri(ilon,ilev) )
     67                 rkoho3 = 1.7e-12*exp( -940./t_seri(ilon,ilev) )
     68                 
     69                 ! rkso2oh= so2 + oh + m  ->hso3 + m-> h2so4 + m
     70                 ! air density M(molecule/cm3) = Pa/(K*1.38e-19)
     71                 rmairdens =  pplay(ilon,ilev)/(t_seri(ilon,ilev)*1.38e-19)
     72                 ! JPL 2020
     73                 rrak0 = 3.3e-31*( (300./t_seri(ilon,ilev))**4.3 )
     74                 rrak1 = 1.6e-12
     75                 rrate = (rrak0*rmairdens) / ( 1. + rrak0*rmairdens/rrak1 )
     76                 rreduce = 1./( 1. + log10((rrak0*rmairdens)/rrak1)**2 )
     77                 rkso2oh = rrate*(0.6**rreduce)
     78                 
     79                 ! O3 (molec/cm3): O3_clim in vmr
     80                 rrak0 = O3_clim(ilon,ilev)*rmairdens
     81                 ! SO2 (molec/cm3): convert from kg/kgA
     82                 rrak1 = tr_seri(ilon,ilev,id_SO2_strat) &
     83                      & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mSO2mol
     84                 
     85                 IF (rrak1 .GE. 0.0) THEN
     86                    rreduce =( (rkho2o3+rkoho3)*rrak0 + rkso2oh*rrak1 ) / &
     87                         (  (rkho2o3+rkoho3)*rrak0 )
     88                    ! reduce OH, so extend SO2 lifetime
     89                    rreduce = rreduce*SO2_lifetime(ilon,ilev)
     90                 ELSE
     91                    rreduce =  SO2_lifetime(ilon,ilev)
     92                 ENDIF
     93                 ! end of IF (rrak1) THEN
     94              ELSE
     95                 rreduce =  SO2_lifetime(ilon,ilev)
     96              ENDIF
     97              ! end of IF (flag_OH_reduced) THEN
     98             
     99              ! Check lifetime rreduce < timestep*1.5 (such as SO2 loss > 0.5*SO2) with exp(-1/1.5)=0.52
     100              ! Check lifetime rreduce < timestep*3 (such as SO2 loss > 0.28*SO2) with exp(-1/3)=0.72
     101              IF(flag_min_rreduce) THEN
     102                 IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys
     103              ENDIF
     104              budg_3D_so2_to_h2so4(ilon,ilev)=tr_seri(ilon,ilev,id_SO2_strat)*(1.0-exp(-pdtphys/rreduce))
     105              tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - budg_3D_so2_to_h2so4(ilon,ilev)
     106              tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) + mH2SO4mol/mSO2mol*budg_3D_so2_to_h2so4(ilon,ilev)
     107           ENDIF
     108           ! IF (SO2_lifetime(ilon,ilev).GT.0.0 .AND. SO2_lifetime(ilon,ilev).LT.1.E10) THEN
     109           
     110           
     111           IF (flag_H2SO4_photolysis) THEN
     112              !--convert H2SO4 to SO2 using HCl photolysis (Rinsland et al, 1995)
     113              ! 4.6e-8 sec-1: J strato taken from Feierabend, et al., Chem. Phys. Lett., 2006.
     114              ! error before:  ...exp(-pdtphys/4.6e-8) so negligible photolysis,
     115              !                 should have been exp(-pdtphys*4.6e-8)
     116              ! Also: Lane and Kjaergaard: ,J. Phys. Chem. A, 2008
     117              ! We find that the dominant photodissociation mechanism of H2SO4
     118              ! below 70 km is absorption in the visible region by OH stretching overtone
     119              ! transitions, whereas above 70 km, absorption of Lyman-R ...
     120              ! In 2D model, 0.3 Rinsland et al, 1995 but Burkholder, Mills and McKeen say 0.016
     121              ! JH2SO4=JHCL*0.3   (*0.016 would be too slow)
     122              ! test: quick sequential integration for SO2 to H2SO4 and reverse
     123             
     124              IF (H2SO4_lifetime(ilon,ilev).GT.0.0 .AND. H2SO4_lifetime(ilon,ilev).LT.1.E10) THEN
     125                 
     126                 rreduce = H2SO4_lifetime(ilon,ilev)
     127                 dummyso4toso2 = 0.0
     128                 
     129                 ! Check lifetime rreduce < timestep*1.5 (such as H2SO4 loss > 0.5*H2SO4) with exp(-1/1.5)=0.52
     130                 ! Check lifetime rreduce < timestep*3 (such as H2SO4 loss > 0.28*H2SO4) with exp(-1/3)=0.72
     131                 IF(flag_min_rreduce) THEN
     132                    IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys
     133                 ENDIF
     134                 dummyso4toso2 = (mSO2mol/mH2SO4mol)*tr_seri(ilon,ilev,id_H2SO4_strat)*(1.0-exp(-pdtphys/rreduce))
     135                 budg_3D_so2_to_h2so4(ilon,ilev) = budg_3D_so2_to_h2so4(ilon,ilev) + dummyso4toso2
     136                 tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + dummyso4toso2
     137                 tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) - (mH2SO4mol/mSO2mol)*dummyso4toso2
     138              ENDIF
     139              ! IF (H2SO4_lifetime(ilon,ilev).GT.0.0 .AND. H2SO4_lifetime(ilon,ilev).LT.1.E10) THEN
     140           ENDIF
     141           ! end of IF (flag_H2SO4_photolysis) THEN
     142           
     143           !convert budget from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic
     144           budg_3D_so2_to_h2so4(ilon,ilev)=budg_3D_so2_to_h2so4(ilon,ilev)*mSatom/mSO2mol* &
     145                (paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
     146           budg_so2_to_h2so4(ilon)=budg_so2_to_h2so4(ilon)+budg_3D_so2_to_h2so4(ilon,ilev)
     147        ENDIF
     148        ! IF (is_strato(ilon,ilev)) THEN
     149     ENDDO ! DO (klev)
     150  ENDDO ! DO (klon)
     151 
    47152END SUBROUTINE SO2_TO_H2SO4

  • LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90

    r4513 r4601  
    2323    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
    2424    USE mod_phys_lmdz_para, only: gather, scatter
    25     USE phys_cal_mod, ONLY : year_len, mth_len, year_cur, mth_cur, day_cur, hour
     25    USE phys_cal_mod, ONLY : year_len, year_cur, mth_cur, day_cur, hour
    2626    USE sulfate_aer_mod
    2727    USE phys_local_var_mod, ONLY: stratomask
    2828    USE YOMCST
    2929    USE print_control_mod, ONLY: lunout
    30     USE strataer_mod
    31 
     30    USE strataer_local_var_mod
     31   
    3232    IMPLICIT NONE
    3333
    … …  
    5858!----------------
    5959    REAL                                   :: m_aer_emiss_vol_daily ! daily injection mass emission
     60    REAL                                   :: m_aer               ! aerosol mass
    6061    INTEGER                                :: it, k, i, ilon, ilev, itime, i_int, ieru
    6162    LOGICAL,DIMENSION(klon,klev)           :: is_strato           ! true = above tropopause, false = below
    … …  
    7879    REAL                                   :: theta_min, theta_max ! for SAI computation between two latitudes
    7980    REAL                                   :: dlat_loc
     81    REAL                                   :: latmin,latmax,lonmin,lonmax ! lat/lon min/max for injection
     82    REAL                                   :: sigma_alt, altemiss ! injection altitude + sigma for distrib
     83    REAL                                   :: pdt,stretchlong     ! physic timestep, stretch emission over one day
     84   
    8085    INTEGER                                :: injdur_sai          ! injection duration for SAI case [days]
    8186    INTEGER                                :: yr, is_bissext
    8287
    83     IF (is_mpi_root) THEN
     88    IF (is_mpi_root .AND. flag_verbose_strataer) THEN
    8489       WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
    85        WRITE(lunout,*) 'IN traccoag flag_sulf_emit: ',flag_sulf_emit
     90       WRITE(lunout,*) 'IN traccoag flag_emit: ',flag_emit
    8691    ENDIF
    8792   
    … …  
    126131!--calculate mass of air in every grid box
    127132    DO ilon=1, klon
    128     DO ilev=1, klev
    129       m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
     133       DO ilev=1, klev
     134          m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
     135       ENDDO
    130136    ENDDO
    131     ENDDO
    132 
    133 !    IF (debutphy) THEN
    134 !      CALL gather(tr_seri, tr_seri_glo)
    135 !      IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN
    136 !--initialising tracer concentrations to zero
    137 !        DO it=1, nbtr
    138 !        tr_seri(:,:,it)=0.0
    139 !        ENDDO
    140 !      ENDIF
    141 !    ENDIF
    142 
     137   
    143138!--initialise emission diagnostics
     139    budg_emi(:,1)=0.0
    144140    budg_emi_ocs(:)=0.0
    145141    budg_emi_so2(:)=0.0
    … …  
    147143    budg_emi_part(:)=0.0
    148144
    149 !--sulfur emission, depending on chosen scenario (flag_sulf_emit)
    150     SELECT CASE(flag_sulf_emit)
     145!--sulfur emission, depending on chosen scenario (flag_emit)
     146    SELECT CASE(flag_emit)
    151147
    152148    CASE(0) ! background aerosol
    … …  
    159155          IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
    160156               day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN
    161              !
    162              ! daily injection mass emission - NL
    163              m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
    164              WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), &
     157
     158             ! daily injection mass emission
     159             m_aer=m_aer_emiss_vol(ieru,1)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
     160             !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
     161             m_aer=m_aer*(mSO2mol/mSatom)
     162             
     163             WRITE(lunout,*) 'IN traccoag m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru,1), &
    165164                  'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', &
    166                   (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily,'ieru=',ieru
     165                  (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer,'ieru=',ieru
    167166             WRITE(lunout,*) 'IN traccoag, dlon=',dlon
    168              DO i=1,klon
    169                 !Pinatubo eruption at 15.14N, 120.35E
    170                 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
    171                 WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc
    172                 IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. &
    173                      xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN
    174                    !
    175                    WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur
    176                    WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily
    177                    !         compute altLMDz
    178                    altLMDz(:)=0.0
    179                    DO k=1, klev
    180                       zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
    181                       zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
    182                       zdz=zdm(k)/zrho                      !thickness of layer in m
    183                       altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
    184                    ENDDO
    185 
    186                    SELECT CASE(flag_sulf_emit_distrib)
    187                    
    188                    CASE(0) ! Gaussian distribution
    189                    !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
    190                    f_lay_sum=0.0
    191                    DO k=1, klev
    192                       f_lay_emiss(k)=0.0
    193                       DO i_int=1, n_int_alt
    194                          alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
    195                          f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol(ieru))* &
    196                               &              exp(-0.5*((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))**2.)*   &
    197                               &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
    198                       ENDDO
    199                       f_lay_sum=f_lay_sum+f_lay_emiss(k)
    200                    ENDDO
    201                    
    202                    CASE(1) ! Uniform distribution
    203                    ! In this case, parameter sigma_alt_vol(ieru) is considered to be half the
    204                    ! height of the injection, centered around altemiss_vol(ieru)
    205                    DO k=1, klev
    206                       f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- &
    207                       & max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru))
    208                       f_lay_sum=f_lay_sum+f_lay_emiss(k)
    209                    ENDDO
    210 
    211                    END SELECT        ! End CASE over flag_sulf_emit_distrib)
    212 
    213                    WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru)
    214                    WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss
    215                    !correct for step integration error
    216                    f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
    217                    !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
    218                    !vertically distributed emission
    219                    DO k=1, klev
    220                       ! stretch emission over one day of Pinatubo eruption
    221                       emission=m_aer_emiss_vol_daily*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys)
    222                       tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
    223                       budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
    224                    ENDDO
    225                 ENDIF ! emission grid cell
    226              ENDDO ! klon loop
    227              WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily
     167             
     168             latmin=xlat_min_vol(ieru)
     169             latmax=xlat_max_vol(ieru)
     170             lonmin=xlon_min_vol(ieru)
     171             lonmax=xlon_max_vol(ieru)
     172             altemiss = altemiss_vol(ieru)
     173             sigma_alt = sigma_alt_vol(ieru)
     174             pdt=pdtphys
     175             ! stretch emission over one day of eruption
     176             stretchlong = 1.
     177             
     178             CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,tr_seri,&
     179                  m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
     180                  stretchlong,1,0)
     181             
    228182          ENDIF ! emission period
    229183       ENDDO ! eruption number
    … …  
    305259        .AND. ( day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur ) &
    306260        .AND. ( day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end )) THEN
    307 
    308      DO i=1,klon
    309         dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
    310         IF  ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. &
    311           &   xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
    312 !
    313 !         compute altLMDz
    314           altLMDz(:)=0.0
    315           DO k=1, klev
    316             zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
    317             zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
    318             zdz=zdm(k)/zrho                      !thickness of layer in m
    319             altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
    320           ENDDO
    321 
    322           SELECT CASE(flag_sulf_emit_distrib)
    323 
    324           CASE(0) ! Gaussian distribution
    325           !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
    326           f_lay_sum=0.0
    327                DO k=1, klev
    328                      f_lay_emiss(k)=0.0
    329                      DO i_int=1, n_int_alt
    330                          alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
    331                          f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
    332                          &              exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   &
    333                          &              (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
    334                      ENDDO
    335                      f_lay_sum=f_lay_sum+f_lay_emiss(k)
    336                ENDDO
    337 
    338           CASE(1) ! Uniform distribution
    339           f_lay_sum=0.0
    340           ! In this case, parameter sigma_alt_vol(ieru) is considered to be half
    341           ! the height of the injection, centered around altemiss_sai
    342                DO k=1, klev
    343                     f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
    344                     & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
    345                     f_lay_sum=f_lay_sum+f_lay_emiss(k)
    346                ENDDO
    347 
    348           END SELECT ! Gaussian or uniform distribution
    349 
    350           !correct for step integration error
    351           f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
    352           !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
    353           !vertically distributed emission
    354           DO k=1, klev
    355             ! stretch emission over whole year (360d)
    356             emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(injdur_sai)/86400. 
    357             tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
    358             budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
    359           ENDDO
    360 
    361 !          !emission as monodisperse particles with 0.1um dry radius (BIN21)
    362 !          !vertically distributed emission
    363 !          DO k=1, klev
    364 !            ! stretch emission over whole year (360d)
    365 !            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(year_len)/86400.
    366 !            tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys
    367 !            budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol
    368 !          ENDDO
    369         ENDIF ! emission grid cell
    370       ENDDO ! klon loop
    371 
     261       
     262       m_aer=m_aer_emiss_sai
     263       !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
     264       m_aer=m_aer*(mSO2mol/mSatom)
     265       
     266       latmin=xlat_sai
     267       latmax=xlat_sai
     268       lonmin=xlon_sai
     269       lonmax=xlon_sai
     270       altemiss = altemiss_sai
     271       sigma_alt = sigma_alt_sai
     272       pdt=0.
     273       ! stretch emission over whole year (360d)
     274       stretchlong=FLOAT(year_len)
     275       
     276       CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
     277            m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
     278            stretchlong,1,0)
     279       
     280       budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
    372281     ENDIF ! Condition over injection dates
    373282
    … …  
    375284            !     lat_min and lat_max
    376285
    377     DO i=1,klon
    378 !       SAI scenario with continuous emission
    379         dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
    380         theta_min = max(xlat(i)-dlat_loc,xlat_min_sai)
    381         theta_max = min(xlat(i)+dlat_loc,xlat_max_sai)
    382         IF  ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. &
    383           &   xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
    384 !
    385 !         compute altLMDz
    386           altLMDz(:)=0.0
    387           DO k=1, klev
    388             zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
    389             zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
    390             zdz=zdm(k)/zrho                      !thickness of layer in m
    391             altLMDz(k+1)=altLMDz(k)+zdz          !altitude of interface
    392           ENDDO
    393 
    394           SELECT CASE(flag_sulf_emit_distrib)
    395 
    396           CASE(0) ! Gaussian distribution
    397           !compute distribution of emission to vertical model layers (based on
    398           !Gaussian peak in altitude)
    399           f_lay_sum=0.0
    400                DO k=1, klev
    401                      f_lay_emiss(k)=0.0
    402                      DO i_int=1, n_int_alt
    403                          alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
    404                          f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* &
    405                          & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)*   &
    406                          & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
    407                      ENDDO
    408                      f_lay_sum=f_lay_sum+f_lay_emiss(k)
    409                ENDDO
    410 
    411           CASE(1) ! Uniform distribution
    412           f_lay_sum=0.0
    413           ! In this case, parameter sigma_alt_vol(ieru) is considered to be half
    414           ! the height of the injection, centered around altemiss_sai
    415                DO k=1, klev
    416                     f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
    417                     & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
    418                     f_lay_sum=f_lay_sum+f_lay_emiss(k)
    419                ENDDO
    420 
    421           END SELECT ! Gaussian or uniform distribution
    422 
    423           !correct for step integration error
    424           f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
    425           !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
    426           !vertically distributed emission
    427           DO k=1, klev
    428             ! stretch emission over whole year (360d)
    429             emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/ &
    430                       & FLOAT(year_len)/86400.*(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ &
    431                       & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI))
    432             tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
    433             budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol
    434           ENDDO
    435 
    436 !          !emission as monodisperse particles with 0.1um dry radius (BIN21)
    437 !          !vertically distributed emission
    438 !          DO k=1, klev
    439 !            ! stretch emission over whole year (360d)
    440 !            emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400
    441 !            tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys
    442 !            budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol
    443 !          ENDDO
    444         ENDIF ! emission grid cell
    445       ENDDO ! klon loop
    446 
    447     END SELECT ! emission scenario (flag_sulf_emit)
     286       m_aer=m_aer_emiss_sai
     287       !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
     288       m_aer=m_aer*(mSO2mol/mSatom)
     289
     290       latmin=xlat_min_sai
     291       latmax=xlat_max_sai
     292       lonmin=xlon_sai
     293       lonmax=xlon_sai
     294       altemiss = altemiss_sai
     295       sigma_alt = sigma_alt_sai
     296       pdt=0.
     297       ! stretch emission over whole year (360d)
     298       stretchlong=FLOAT(year_len)
     299
     300       CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
     301            m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
     302            stretchlong,1,0)
     303
     304       budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
     305       
     306    END SELECT ! emission scenario (flag_emit)
    448307
    449308!--read background concentrations of OCS and SO2 and lifetimes from input file
    … …  
    463322    CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
    464323
    465 !--call sedimentation routine 
     324!--call sedimentation routine
    466325    CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
    467326
    … …  
    480339      ENDDO
    481340    ENDDO
    482 
    483 !    CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2')
    484 !    DO it=1, nbtr_bin
    485 !      CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ')
    486 !    ENDDO
    487 
     341   
    488342  END SUBROUTINE traccoag
    489343

  • LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90

    r4575 r4601  
    542542!
    543543! variables for stratospheric aerosol
     544      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: d_q_emiss
     545!$OMP THREADPRIVATE(d_q_emiss)
    544546      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: R2SO4
    545547!$OMP THREADPRIVATE(R2SO4)
    … …  
    556558      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: SO2_lifetime
    557559!$OMP THREADPRIVATE(SO2_lifetime)
     560      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: H2SO4_lifetime
     561!$OMP THREADPRIVATE(H2SO4_lifetime)
     562      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: O3_clim
     563!$OMP THREADPRIVATE(O3_clim)
    558564      REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: alpha_bin
    559565!$OMP THREADPRIVATE(alpha_bin)
    … …  
    920926
    921927#ifdef CPP_StratAer
     928      ALLOCATE (d_q_emiss(klon,klev))
    922929      ALLOCATE (R2SO4(klon,klev))
    923930      ALLOCATE (DENSO4(klon,klev))
    … …  
    933940      ALLOCATE (OCS_lifetime(klon,klev))
    934941      ALLOCATE (SO2_lifetime(klon,klev))
     942      ALLOCATE (H2SO4_lifetime(klon,klev))
     943      ALLOCATE (O3_clim(klon,klev))
    935944      ALLOCATE (alpha_bin(nbands_sw_rrtm+nbands_lw_rrtm+nwave,nbtr))
    936945      ALLOCATE (piz_bin(nbands_sw_rrtm+nbands_lw_rrtm+nwave,nbtr))
    … …  
    12271236#ifdef CPP_StratAer
    12281237! variables for strat. aerosol CK
     1238      DEALLOCATE (d_q_emiss)
    12291239      DEALLOCATE (R2SO4)
    12301240      DEALLOCATE (DENSO4)
    … …  
    12341244      DEALLOCATE (SO2_lifetime)
    12351245      DEALLOCATE (OCS_lifetime)
     1246      DEALLOCATE (H2SO4_lifetime)
     1247      DEALLOCATE (O3_clim)
    12361248      DEALLOCATE (alpha_bin)
    12371249      DEALLOCATE (piz_bin)

  • LMDZ6/trunk/libf/phylmd/physiq_mod.F90

    r4596 r4601  
    126126
    127127#ifdef CPP_StratAer
    128     USE strataer_mod, ONLY: strataer_init
     128    USE phys_local_var_mod, ONLY: d_q_emiss
     129    USE strataer_local_var_mod
     130    USE strataer_nuc_mod, ONLY: strataer_nuc_init
     131    USE strataer_emiss_mod, ONLY: strataer_emiss_init
    129132#endif
    130133
    … …  
    13481351#ifdef CPP_StratAer
    13491352       CALL strataer_init
     1353       CALL strataer_nuc_init
     1354       CALL strataer_emiss_init
    13501355#endif
    13511356
    … …  
    15311536          piz_aero(:,:,:,:) = 0.
    15321537          cg_aero(:,:,:,:) = 0.
     1538          d_q_ch4(:,:) = 0.
    15331539
    15341540          config_inca='none' ! default
    … …  
    48524858    !
    48534859    !
     4860#ifdef CPP_StratAer
     4861    IF (ok_qemiss) THEN
     4862       flh2o=1
     4863       IF(flag_verbose_strataer) THEN
     4864          print *,'IN physiq_mod: ok_qemiss =yes (',ok_qemiss,'), flh2o=',flh2o
     4865          print *,'IN physiq_mod: flag_emit=',flag_emit,', nErupt=',nErupt
     4866          print *,'IN physiq_mod: nAerErupt=',nAerErupt
     4867       ENDIF
     4868       
     4869       SELECT CASE(flag_emit)
     4870       CASE(1) ! emission volc H2O dans LMDZ
     4871          DO ieru=1, nErupt
     4872             IF (year_cur==year_emit_vol(ieru).AND.&
     4873                  mth_cur==mth_emit_vol(ieru).AND.&
     4874                  day_cur>=day_emit_vol(ieru).AND.&
     4875                  day_cur<(day_emit_vol(ieru)+injdur)) THEN
     4876               
     4877                IF(flag_verbose_strataer) print *,'IN physiq_mod: date=',year_cur,mth_cur,day_cur
     4878                ! initialisation tendance q emission
     4879                d_q_emiss(:,:)=0.
     4880                ! daily injection mass emission - NL
     4881                m_H2O_emiss_vol_daily = m_H2O_emiss_vol(ieru)/(REAL(injdur)&
     4882                     *REAL(ponde_lonlat_vol(ieru)))
     4883                !
     4884                CALL STRATEMIT(pdtphys,pdtphys,latitude_deg,longitude_deg,t_seri,&
     4885                    pplay,paprs,tr_seri,&
     4886                    m_H2O_emiss_vol_daily,&
     4887                    xlat_min_vol(ieru),xlat_max_vol(ieru),&
     4888                    xlon_min_vol(ieru),xlon_max_vol(ieru),&
     4889                    altemiss_vol(ieru),sigma_alt_vol(ieru),1,1.,&
     4890                    nAerErupt+1,0)
     4891               
     4892                IF(flag_verbose_strataer) print *,'IN physiq_mod: min max d_q_emiss=',&
     4893                     minval(d_q_emiss),maxval(d_q_emiss)
     4894               
     4895                CALL add_phys_tend(du0, dv0, dt0, d_q_emiss, dql0, dqi0, dqbs0, paprs, &
     4896                     'q_emiss',abortphy,flag_inhib_tend,itap,0)
     4897                IF (abortphy==1) Print*,'ERROR ABORT TEND EMISS'
     4898             ENDIF
     4899          ENDDO
     4900          flh2o=0
     4901       END SELECT ! emission scenario (flag_emit)
     4902    ENDIF
     4903#endif
    48544904
    48554905!===============================================================
    … …  
    53315381      qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
    53325382    ENDDO
     5383
     5384#ifdef CPP_StratAer
     5385    IF (ok_qemiss) THEN
     5386       DO k = 1, klev
     5387          qql1(:) = qql1(:)+d_q_emiss(:,k)*zmasse(:,k)
     5388       ENDDO
     5389    ENDIF
     5390#endif
     5391    IF (ok_qch4) THEN
     5392       DO k = 1, klev
     5393          qql1(:) = qql1(:)+d_q_ch4_dtime(:,k)*zmasse(:,k)
     5394       ENDDO
     5395    ENDIF
     5396   
    53335397    DO i = 1, klon
    53345398      !--compute ratio of what q+ql should be with conservation to what it is

  • LMDZ6/trunk/libf/phylmd/phytrac_mod.F90

    r4590 r4601  
    147147    USE phys_local_var_mod, ONLY: budg_dep_dry_part,  budg_dep_wet_part
    148148    USE infotrac_phy, ONLY: nbtr_sulgas, id_OCS_strat, id_SO2_strat, id_H2SO4_strat
     149    USE strataer_nuc_mod, ONLY : tracstrataer_init
    149150    USE aerophys
    150151#endif
    … …  
    515516       ELSE IF (type_trac == 'coag') THEN
    516517          source(:,:)=0.
    517           DO it= 1, nbtr_sulgas
    518             aerosol(it)=.FALSE.
    519             IF (it==id_H2SO4_strat) aerosol(it)=.TRUE.
    520           ENDDO
    521           DO it= nbtr_sulgas+1, nbtr
    522             aerosol(it)=.TRUE.
    523           ENDDO
     518          CALL tracstrataer_init(aerosol,lessivage)       ! init aerosols and lessivage param
    524519#endif
    525520       ELSE IF (type_trac == 'lmdz') THEN

  • LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90

    r4594 r4601  
    126126
    127127#ifdef CPP_StratAer
    128     USE strataer_mod, ONLY: strataer_init
     128    USE strataer_local_var_mod
     129    USE strataer_nuc_mod, ONLY: strataer_nuc_init
     130    USE strataer_emiss_mod, ONLY: strataer_emiss_init
    129131#endif
    130132
    … …  
    14141416#ifdef CPP_StratAer
    14151417       CALL strataer_init
     1418       CALL strataer_nuc_init
     1419       CALL strataer_emiss_init
    14161420#endif
    14171421
    … …  
    15991603          piz_aero(:,:,:,:) = 0.
    16001604          cg_aero(:,:,:,:) = 0.
     1605          d_q_ch4(:,:) = 0.
    16011606
    16021607          config_inca='none' ! default
    … …  
    16101615       piz_aero(:,:,:,:) = 1.
    16111616       cg_aero(:,:,:,:)  = 0.
     1617       d_q_ch4(:,:) = 0.
    16121618
    16131619       IF (aerosol_couple .AND. (config_inca /= "aero" &
    … …  
    61546160    !
    61556161    !
     6162#ifdef CPP_StratAer
     6163    IF (ok_qemiss) THEN
     6164       flh2o=1
     6165       IF(flag_verbose_strataer) THEN
     6166          print *,'IN physiq_mod: ok_qemiss =yes (',ok_qemiss,'), flh2o=',flh2o
     6167          print *,'IN physiq_mod: flag_emit=',flag_emit,', nErupt=',nErupt
     6168          print *,'IN physiq_mod: nAerErupt=',nAerErupt
     6169       ENDIF
     6170       
     6171       SELECT CASE(flag_emit)
     6172       CASE(1) ! emission volc H2O dans LMDZ
     6173          DO ieru=1, nErupt
     6174             IF (year_cur==year_emit_vol(ieru).AND.&
     6175                  mth_cur==mth_emit_vol(ieru).AND.&
     6176                  day_cur>=day_emit_vol(ieru).AND.&
     6177                  day_cur<(day_emit_vol(ieru)+injdur)) THEN
     6178               
     6179                IF(flag_verbose_strataer) print *,'IN physiq_mod: date=',year_cur,mth_cur,day_cur
     6180                ! initialisation tendance q emission
     6181                d_q_emiss(:,:)=0.
     6182                ! daily injection mass emission - NL
     6183                m_H2O_emiss_vol_daily = m_H2O_emiss_vol(ieru)/(REAL(injdur)&
     6184                     *REAL(ponde_lonlat_vol(ieru)))
     6185                !
     6186                CALL STRATEMIT(pdtphys,pdtphys,latitude_deg,longitude_deg,t_seri,&
     6187                     pplay,paprs,tr_seri,&
     6188                     m_H2O_emiss_vol_daily,&
     6189                     xlat_min_vol(ieru),xlat_max_vol(ieru),&
     6190                     xlon_min_vol(ieru),xlon_max_vol(ieru),&
     6191                     altemiss_vol(ieru),&
     6192                     sigma_alt_vol(ieru),1,&
     6193                     1,nAerErupt+1,0)
     6194               
     6195                IF(flag_verbose_strataer) print *,'IN physiq_mod: min max d_q_emiss=',&
     6196                     minval(d_q_emiss),maxval(d_q_emiss)
     6197               
     6198                CALL add_phys_tend(du0, dv0, dt0, d_q_emiss, dql0, dqi0, dqbs0, paprs, &
     6199                     'q_emiss',abortphy,flag_inhib_tend,itap,0)
     6200                IF (abortphy==1) Print*,'ERROR ABORT TEND EMISS'
     6201             ENDIF
     6202          ENDDO
     6203          flh2o=0
     6204       END SELECT ! emission scenario (flag_emit)
     6205    ENDIF
     6206#endif
    61566207
    61576208!===============================================================
    … …  
    67446795      qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
    67456796    ENDDO
     6797   
     6798#ifdef CPP_StratAer
     6799    IF (ok_qemiss) THEN
     6800       DO k = 1, klev
     6801          qql1(:) = qql1(:)+d_q_emiss(:,k)*zmasse(:,k)
     6802       ENDDO
     6803    ENDIF
     6804#endif
     6805    IF (ok_qch4) THEN
     6806       DO k = 1, klev
     6807          qql1(:) = qql1(:)+d_q_ch4_dtime(:,k)*zmasse(:,k)
     6808       ENDDO
     6809    ENDIF
     6810   
    67466811    DO i = 1, klon
    67476812      !--compute ratio of what q+ql should be with conservation to what it is
Note: See TracChangeset for help on using the changeset viewer.