Ignore:
Timestamp:
Jul 24, 2024, 12:17:33 PM (7 months ago)
Author:
abarral
Message:

Put abort_physic into a module
Remove -g option from makelmdz_fcm, since that option is linked to a header file that isn't included anywhere.
(lint) light lint on traversed files

Location:
LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/strataer_emiss_mod.F90

    r5110 r5111  
    1010    USE print_control_mod, ONLY: lunout
    1111    USE lmdz_phys_para, ONLY: is_master
     12    USE lmdz_abort_physic, ONLY: abort_physic
    1213
    1314    ! Local variables
    14     INTEGER                  :: ispec
    15 
    16     WRITE(lunout,*) 'IN STRATAER_EMISS INIT WELCOME!'
     15    INTEGER :: ispec
     16
     17    WRITE(lunout, *) 'IN STRATAER_EMISS INIT WELCOME!'
    1718
    1819    IF (flag_emit==1 .OR. flag_emit==4) THEN ! Volcano
    19        CALL getin_p('nErupt',nErupt) !eruption nb
    20        CALL getin_p('injdur',injdur) !injection duration
    21 
    22        IF (flag_emit==1) THEN
    23           CALL getin_p('nAerErupt', nAerErupt) !sulfur aer nb
    24        ELSEIF (flag_emit==4) THEN
    25           CALL getin_p('nSpeciesErupt', nSpeciesErupt) !chimical species nb
    26        ENDIF
    27 
    28        IF (nErupt>0) THEN
    29           ALLOCATE(year_emit_vol(nErupt),mth_emit_vol(nErupt),day_emit_vol(nErupt))
    30           year_emit_vol=0 ; mth_emit_vol=0 ; day_emit_vol=0
    31           ALLOCATE(altemiss_vol(nErupt),sigma_alt_vol(nErupt))
    32           ALLOCATE(xlat_min_vol(nErupt),xlon_min_vol(nErupt))
    33           ALLOCATE(xlat_max_vol(nErupt),xlon_max_vol(nErupt))
    34           altemiss_vol=0. ; sigma_alt_vol=0.
    35           xlon_min_vol=0. ; xlon_max_vol=0.
    36           xlat_min_vol=0. ; xlat_max_vol=0.
    37           ! injection params (dates, loc, injections params)
    38           CALL getin_p('year_emit_vol',year_emit_vol)
    39           CALL getin_p('mth_emit_vol',mth_emit_vol)
    40           CALL getin_p('day_emit_vol',day_emit_vol)
    41           CALL getin_p('altemiss_vol',altemiss_vol)
    42           CALL getin_p('sigma_alt_vol',sigma_alt_vol)
    43           CALL getin_p('xlon_min_vol',xlon_min_vol)
    44           CALL getin_p('xlon_max_vol',xlon_max_vol)
    45           CALL getin_p('xlat_min_vol',xlat_min_vol)
    46           CALL getin_p('xlat_max_vol',xlat_max_vol)
    47           IF (flag_emit==1) THEN
    48              ALLOCATE(m_sulf_emiss_vol(nErupt))
    49              ALLOCATE(m_aer_emiss_vol(nErupt,nAerErupt))
    50              m_aer_emiss_vol=0. ; m_sulf_emiss_vol=0.
    51              IF (ok_qemiss) then
    52                 ALLOCATE(m_H2O_emiss_vol(nErupt))
    53                 ALLOCATE(m_H2O_emiss_vol_daily(nErupt))
    54                 !                ALLOCATE(d_q_emiss(klon,klev))
    55                 ALLOCATE(budg_emi(klon,nAerErupt+1))
    56                 m_H2O_emiss_vol(:)=0.
    57                 m_H2O_emiss_vol_daily(:)=0.
    58                 !                d_q_emiss(:,:)=0.
    59              ELSE
    60                 ALLOCATE(budg_emi(klon,nAerErupt))
    61              ENDIF
    62              budg_emi(:,:)=0.
    63           ELSEIF (flag_emit==4) THEN
    64              ALLOCATE(m_Chlore_emiss_vol(nErupt))
    65              ALLOCATE(id_HCl)
    66              ALLOCATE(m_Brome_emiss_vol(nErupt))
    67              ALLOCATE(id_HBr)
    68              ALLOCATE(id_species(nSpeciesErupt))
    69              id_species = 0
    70              ALLOCATE(m_species_emiss_vol(nErupt,nSpeciesErupt))
    71              m_species_emiss_vol=0.
    72              ALLOCATE(m_NOx_emiss_vol(nErupt))
    73              ALLOCATE(m_H2O_emiss_vol(nErupt))
    74              m_Chlore_emiss_vol=0. ; m_Brome_emiss_vol=0.
    75              m_NOx_emiss_vol=0. ; m_H2O_emiss_vol=0.
    76              ALLOCATE(id_NOx)
    77              ALLOCATE(id_H2O)
    78              ALLOCATE(budg_emi(klon,nSpeciesErupt))
    79              budg_emi(:,:)=0.
     20      CALL getin_p('nErupt', nErupt) !eruption nb
     21      CALL getin_p('injdur', injdur) !injection duration
     22
     23      IF (flag_emit==1) THEN
     24        CALL getin_p('nAerErupt', nAerErupt) !sulfur aer nb
     25      ELSEIF (flag_emit==4) THEN
     26        CALL getin_p('nSpeciesErupt', nSpeciesErupt) !chimical species nb
     27      ENDIF
     28
     29      IF (nErupt>0) THEN
     30        ALLOCATE(year_emit_vol(nErupt), mth_emit_vol(nErupt), day_emit_vol(nErupt))
     31        year_emit_vol = 0 ; mth_emit_vol = 0 ; day_emit_vol = 0
     32        ALLOCATE(altemiss_vol(nErupt), sigma_alt_vol(nErupt))
     33        ALLOCATE(xlat_min_vol(nErupt), xlon_min_vol(nErupt))
     34        ALLOCATE(xlat_max_vol(nErupt), xlon_max_vol(nErupt))
     35        altemiss_vol = 0. ; sigma_alt_vol = 0.
     36        xlon_min_vol = 0. ; xlon_max_vol = 0.
     37        xlat_min_vol = 0. ; xlat_max_vol = 0.
     38        ! injection params (dates, loc, injections params)
     39        CALL getin_p('year_emit_vol', year_emit_vol)
     40        CALL getin_p('mth_emit_vol', mth_emit_vol)
     41        CALL getin_p('day_emit_vol', day_emit_vol)
     42        CALL getin_p('altemiss_vol', altemiss_vol)
     43        CALL getin_p('sigma_alt_vol', sigma_alt_vol)
     44        CALL getin_p('xlon_min_vol', xlon_min_vol)
     45        CALL getin_p('xlon_max_vol', xlon_max_vol)
     46        CALL getin_p('xlat_min_vol', xlat_min_vol)
     47        CALL getin_p('xlat_max_vol', xlat_max_vol)
     48        IF (flag_emit==1) THEN
     49          ALLOCATE(m_sulf_emiss_vol(nErupt))
     50          ALLOCATE(m_aer_emiss_vol(nErupt, nAerErupt))
     51          m_aer_emiss_vol = 0. ; m_sulf_emiss_vol = 0.
     52          IF (ok_qemiss) then
     53            ALLOCATE(m_H2O_emiss_vol(nErupt))
     54            ALLOCATE(m_H2O_emiss_vol_daily(nErupt))
     55            !                ALLOCATE(d_q_emiss(klon,klev))
     56            ALLOCATE(budg_emi(klon, nAerErupt + 1))
     57            m_H2O_emiss_vol(:) = 0.
     58            m_H2O_emiss_vol_daily(:) = 0.
     59            !                d_q_emiss(:,:)=0.
     60          ELSE
     61            ALLOCATE(budg_emi(klon, nAerErupt))
    8062          ENDIF
    81        ELSE
    82           WRITE(lunout,*) 'ERROR : Using flag_emit=1 or 4 (ie Volcanic eruption) but nErupt (',nErupt,') <=0 !'
    83           CALL abort_physic('strataer_emiss_mod', &
    84                'No eruption define in physiq_def (nErupt=0). Add at one eruption or use background condition.',1)
    85        ENDIF ! fin if nerupt
    86 
    87        IF (flag_emit==1) THEN
    88           CALL getin_p('m_sulf_emiss_vol',m_sulf_emiss_vol)
    89           if (ok_qemiss) then
    90              CALL getin_p('m_H2O_emiss_vol',m_H2O_emiss_vol)
    91           endif
    92        ELSEIF (flag_emit==4) THEN
    93           CALL getin_p('id_species',id_species)
    94           CALL getin_p('m_Chlore_emiss_vol',m_Chlore_emiss_vol)
    95           CALL getin_p('id_HCl',id_HCl)
    96           CALL getin_p('m_Brome_emiss_vol',m_Brome_emiss_vol)
    97           CALL getin_p('id_HBr',id_HBr)
    98           CALL getin_p('m_NOx_emiss_vol',m_NOx_emiss_vol)
    99           CALL getin_p('id_NOx',id_NOx)
    100           CALL getin_p('m_H2O_emiss_vol',m_H2O_emiss_vol)
    101           CALL getin_p('id_H2O',id_H2O)
    102        ENDIF
     63          budg_emi(:, :) = 0.
     64        ELSEIF (flag_emit==4) THEN
     65          ALLOCATE(m_Chlore_emiss_vol(nErupt))
     66          ALLOCATE(id_HCl)
     67          ALLOCATE(m_Brome_emiss_vol(nErupt))
     68          ALLOCATE(id_HBr)
     69          ALLOCATE(id_species(nSpeciesErupt))
     70          id_species = 0
     71          ALLOCATE(m_species_emiss_vol(nErupt, nSpeciesErupt))
     72          m_species_emiss_vol = 0.
     73          ALLOCATE(m_NOx_emiss_vol(nErupt))
     74          ALLOCATE(m_H2O_emiss_vol(nErupt))
     75          m_Chlore_emiss_vol = 0. ; m_Brome_emiss_vol = 0.
     76          m_NOx_emiss_vol = 0. ; m_H2O_emiss_vol = 0.
     77          ALLOCATE(id_NOx)
     78          ALLOCATE(id_H2O)
     79          ALLOCATE(budg_emi(klon, nSpeciesErupt))
     80          budg_emi(:, :) = 0.
     81        ENDIF
     82      ELSE
     83        WRITE(lunout, *) 'ERROR : Using flag_emit=1 or 4 (ie Volcanic eruption) but nErupt (', nErupt, ') <=0 !'
     84        CALL abort_physic('strataer_emiss_mod', &
     85                'No eruption define in physiq_def (nErupt=0). Add at one eruption or use background condition.', 1)
     86      ENDIF ! fin if nerupt
     87
     88      IF (flag_emit==1) THEN
     89        CALL getin_p('m_sulf_emiss_vol', m_sulf_emiss_vol)
     90        if (ok_qemiss) then
     91          CALL getin_p('m_H2O_emiss_vol', m_H2O_emiss_vol)
     92        endif
     93      ELSEIF (flag_emit==4) THEN
     94        CALL getin_p('id_species', id_species)
     95        CALL getin_p('m_Chlore_emiss_vol', m_Chlore_emiss_vol)
     96        CALL getin_p('id_HCl', id_HCl)
     97        CALL getin_p('m_Brome_emiss_vol', m_Brome_emiss_vol)
     98        CALL getin_p('id_HBr', id_HBr)
     99        CALL getin_p('m_NOx_emiss_vol', m_NOx_emiss_vol)
     100        CALL getin_p('id_NOx', id_NOx)
     101        CALL getin_p('m_H2O_emiss_vol', m_H2O_emiss_vol)
     102        CALL getin_p('id_H2O', id_H2O)
     103      ENDIF
    103104
    104105    ELSEIF (flag_emit == 2) THEN ! SAI
    105        CALL getin_p('m_aer_emiss_sai',m_aer_emiss_sai)
    106        CALL getin_p('altemiss_sai',altemiss_sai)
    107        CALL getin_p('sigma_alt_sai',sigma_alt_sai)
    108        CALL getin_p('xlat_sai',xlat_sai)
    109        CALL getin_p('xlon_sai',xlon_sai)
    110        CALL getin_p('year_emit_sai_start',year_emit_sai_start)
    111        CALL getin_p('year_emit_sai_end',year_emit_sai_end)
    112        CALL getin_p('mth_emit_sai_start',mth_emit_sai_start)
    113        CALL getin_p('mth_emit_sai_end',mth_emit_sai_end)
    114        CALL getin_p('day_emit_sai_start',day_emit_sai_start)
    115        CALL getin_p('day_emit_sai_end',day_emit_sai_end)
     106      CALL getin_p('m_aer_emiss_sai', m_aer_emiss_sai)
     107      CALL getin_p('altemiss_sai', altemiss_sai)
     108      CALL getin_p('sigma_alt_sai', sigma_alt_sai)
     109      CALL getin_p('xlat_sai', xlat_sai)
     110      CALL getin_p('xlon_sai', xlon_sai)
     111      CALL getin_p('year_emit_sai_start', year_emit_sai_start)
     112      CALL getin_p('year_emit_sai_end', year_emit_sai_end)
     113      CALL getin_p('mth_emit_sai_start', mth_emit_sai_start)
     114      CALL getin_p('mth_emit_sai_end', mth_emit_sai_end)
     115      CALL getin_p('day_emit_sai_start', day_emit_sai_start)
     116      CALL getin_p('day_emit_sai_end', day_emit_sai_end)
    116117
    117118    ELSEIF (flag_emit == 3) THEN ! SAI between latitudes
    118        CALL getin_p('m_aer_emiss_sai',m_aer_emiss_sai)
    119        CALL getin_p('altemiss_sai',altemiss_sai)
    120        CALL getin_p('sigma_alt_sai',sigma_alt_sai)
    121        CALL getin_p('xlon_sai',xlon_sai)
    122        CALL getin_p('xlat_max_sai',xlat_max_sai)
    123        CALL getin_p('xlat_min_sai',xlat_min_sai)
     119      CALL getin_p('m_aer_emiss_sai', m_aer_emiss_sai)
     120      CALL getin_p('altemiss_sai', altemiss_sai)
     121      CALL getin_p('sigma_alt_sai', sigma_alt_sai)
     122      CALL getin_p('xlon_sai', xlon_sai)
     123      CALL getin_p('xlat_max_sai', xlat_max_sai)
     124      CALL getin_p('xlat_min_sai', xlat_min_sai)
    124125    ENDIF
    125126
    126127    IF (flag_emit == 1) THEN
    127        DO ispec=1,nAerErupt
    128           m_aer_emiss_vol(:,ispec) = m_sulf_emiss_vol(:)
    129        ENDDO
     128      DO ispec = 1, nAerErupt
     129        m_aer_emiss_vol(:, ispec) = m_sulf_emiss_vol(:)
     130      ENDDO
    130131    ELSEIF (flag_emit== 4) THEN
    131        DO ispec=1,nSpeciesErupt
    132           IF(id_species(ispec) == id_HCl) THEN
    133              m_species_emiss_vol(:,ispec) = m_Chlore_emiss_vol(:)
    134           ENDIF
    135           IF (id_species(ispec) == id_HBr) THEN
    136              m_species_emiss_vol(:,ispec) = m_Brome_emiss_vol(:)
    137           ENDIF
    138           IF (id_species(ispec) == id_NOx) THEN
    139              m_species_emiss_vol(:,ispec) = m_NOx_emiss_vol(:)
    140           ENDIF
    141           IF (id_species(ispec) == id_H2O) THEN
    142              m_species_emiss_vol(:,ispec) = m_H2O_emiss_vol(:)
    143           ENDIF
    144        ENDDO
     132      DO ispec = 1, nSpeciesErupt
     133        IF(id_species(ispec) == id_HCl) THEN
     134          m_species_emiss_vol(:, ispec) = m_Chlore_emiss_vol(:)
     135        ENDIF
     136        IF (id_species(ispec) == id_HBr) THEN
     137          m_species_emiss_vol(:, ispec) = m_Brome_emiss_vol(:)
     138        ENDIF
     139        IF (id_species(ispec) == id_NOx) THEN
     140          m_species_emiss_vol(:, ispec) = m_NOx_emiss_vol(:)
     141        ENDIF
     142        IF (id_species(ispec) == id_H2O) THEN
     143          m_species_emiss_vol(:, ispec) = m_H2O_emiss_vol(:)
     144        ENDIF
     145      ENDDO
    145146    ENDIF
    146147
    147148    !============= Injection ponderation =============
    148149    IF (flag_emit > 0) THEN
    149        CALL strataer_ponde_init
    150        WRITE(lunout,*) 'IN STRATAER INIT : ponde_lonlat_vol',ponde_lonlat_vol
     150      CALL strataer_ponde_init
     151      WRITE(lunout, *) 'IN STRATAER INIT : ponde_lonlat_vol', ponde_lonlat_vol
    151152    ENDIF
    152153
    153154    !============= Print params =============
    154155    IF (is_master) THEN
    155        IF (nErupt > 0) then
    156           if (flag_emit == 1 .OR. flag_emit == 4) THEN
    157              WRITE(lunout,*) 'IN STRATAER nErupt: ',nErupt
    158              WRITE(lunout,*) 'IN STRATAER injdur: ',injdur
    159              WRITE(lunout,*) 'IN STRATAER nAerErupt: ',nAerErupt
    160 
    161              WRITE(lunout,*) 'IN STRATAER : year_emit_vol',year_emit_vol
    162              WRITE(lunout,*) 'IN STRATAER : mth_emit_vol',mth_emit_vol
    163              WRITE(lunout,*) 'IN STRATAER : day_emit_vol',day_emit_vol
    164              WRITE(lunout,*) 'IN STRATAER : altemiss_vol',altemiss_vol
    165              WRITE(lunout,*) 'IN STRATAER : sigma_alt_vol',sigma_alt_vol
    166              WRITE(lunout,*) 'IN STRATAER : xlat_min_vol',xlat_min_vol
    167              WRITE(lunout,*) 'IN STRATAER : xlat_max_vol',xlat_max_vol
    168              WRITE(lunout,*) 'IN STRATAER : xlon_min_vol',xlon_min_vol
    169              WRITE(lunout,*) 'IN STRATAER : xlon_max_vol',xlon_max_vol
    170              IF (flag_emit==1) THEN
    171                 WRITE(lunout,*) 'IN STRATAEREMISS : m_sulf_emiss_vol',m_sulf_emiss_vol
    172                 WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_vol',m_aer_emiss_vol
    173                 IF (ok_qemiss) then
    174                    WRITE(lunout,*) 'IN STRATAEREMISS : m_H2O_emiss_vol',m_H2O_emiss_vol
    175                 ENDIF
    176              ENDIF
    177           ELSEIF (flag_emit == 2) THEN
    178              WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_sai',m_aer_emiss_sai
    179              WRITE(lunout,*) 'IN STRATAER : altemiss_sai',altemiss_sai
    180              WRITE(lunout,*) 'IN STRATAER : sigma_alt_sai',sigma_alt_sai
    181              WRITE(lunout,*) 'IN STRATAER : xlat_sai',xlat_sai
    182              WRITE(lunout,*) 'IN STRATAER : xlon_sai',xlon_sai
    183           ELSEIF (flag_emit == 3) THEN
    184              WRITE(lunout,*) 'IN STRATAER : m_aer_emiss_sai',m_aer_emiss_sai
    185              WRITE(lunout,*) 'IN STRATAER : altemiss_sai',altemiss_sai
    186              WRITE(lunout,*) 'IN STRATAER : sigma_alt_sai',sigma_alt_sai
    187              WRITE(lunout,*) 'IN STRATAER : year_emit_sai start/end',year_emit_sai_start,year_emit_sai_end
    188              WRITE(lunout,*) 'IN STRATAER : mth_emit_sai start/end',mth_emit_sai_start,mth_emit_sai_end
    189              WRITE(lunout,*) 'IN STRATAER : day_emit_sai start/end',day_emit_sai_start,day_emit_sai_end
    190              WRITE(lunout,*) 'IN STRATAER : xlat_min_sai',xlat_min_sai
    191              WRITE(lunout,*) 'IN STRATAER : xlat_max_sai',xlat_max_sai
    192              WRITE(lunout,*) 'IN STRATAER : xlon_sai',xlon_sai
     156      IF (nErupt > 0) then
     157        if (flag_emit == 1 .OR. flag_emit == 4) THEN
     158          WRITE(lunout, *) 'IN STRATAER nErupt: ', nErupt
     159          WRITE(lunout, *) 'IN STRATAER injdur: ', injdur
     160          WRITE(lunout, *) 'IN STRATAER nAerErupt: ', nAerErupt
     161
     162          WRITE(lunout, *) 'IN STRATAER : year_emit_vol', year_emit_vol
     163          WRITE(lunout, *) 'IN STRATAER : mth_emit_vol', mth_emit_vol
     164          WRITE(lunout, *) 'IN STRATAER : day_emit_vol', day_emit_vol
     165          WRITE(lunout, *) 'IN STRATAER : altemiss_vol', altemiss_vol
     166          WRITE(lunout, *) 'IN STRATAER : sigma_alt_vol', sigma_alt_vol
     167          WRITE(lunout, *) 'IN STRATAER : xlat_min_vol', xlat_min_vol
     168          WRITE(lunout, *) 'IN STRATAER : xlat_max_vol', xlat_max_vol
     169          WRITE(lunout, *) 'IN STRATAER : xlon_min_vol', xlon_min_vol
     170          WRITE(lunout, *) 'IN STRATAER : xlon_max_vol', xlon_max_vol
     171          IF (flag_emit==1) THEN
     172            WRITE(lunout, *) 'IN STRATAEREMISS : m_sulf_emiss_vol', m_sulf_emiss_vol
     173            WRITE(lunout, *) 'IN STRATAER : m_aer_emiss_vol', m_aer_emiss_vol
     174            IF (ok_qemiss) then
     175              WRITE(lunout, *) 'IN STRATAEREMISS : m_H2O_emiss_vol', m_H2O_emiss_vol
     176            ENDIF
    193177          ENDIF
    194           IF(flag_emit == 4) THEN
    195              WRITE(lunout,*) 'IN STRATAER : nSpeciesErupt: ',nSpeciesErupt
    196              WRITE(lunout,*) 'IN STRATAER : id_species = ',id_species
    197              WRITE(lunout,*) 'IN STRATAER : id_HCl = ',id_HCl
    198              WRITE(lunout,*) 'IN STRATAER : id_HBr = ',id_HBr
    199              WRITE(lunout,*) 'IN STRATAER : id_NOx = ',id_NOx
    200              WRITE(lunout,*) 'IN STRATAER : id_H2O = ',id_H2O
    201              WRITE(lunout,*) 'IN STRATAER : m_Chlore_emiss_vol = ',m_Chlore_emiss_vol
    202              WRITE(lunout,*) 'IN STRATAER : m_Brome_emiss_vol = ',m_Brome_emiss_vol
    203              WRITE(lunout,*) 'IN STRATAER : m_NOx_emiss_vol = ',m_NOx_emiss_vol
    204              WRITE(lunout,*) 'IN STRATAER : m_H2O_emiss_vol = ',m_H2O_emiss_vol
    205           ENDIF
    206        endif
     178        ELSEIF (flag_emit == 2) THEN
     179          WRITE(lunout, *) 'IN STRATAER : m_aer_emiss_sai', m_aer_emiss_sai
     180          WRITE(lunout, *) 'IN STRATAER : altemiss_sai', altemiss_sai
     181          WRITE(lunout, *) 'IN STRATAER : sigma_alt_sai', sigma_alt_sai
     182          WRITE(lunout, *) 'IN STRATAER : xlat_sai', xlat_sai
     183          WRITE(lunout, *) 'IN STRATAER : xlon_sai', xlon_sai
     184        ELSEIF (flag_emit == 3) THEN
     185          WRITE(lunout, *) 'IN STRATAER : m_aer_emiss_sai', m_aer_emiss_sai
     186          WRITE(lunout, *) 'IN STRATAER : altemiss_sai', altemiss_sai
     187          WRITE(lunout, *) 'IN STRATAER : sigma_alt_sai', sigma_alt_sai
     188          WRITE(lunout, *) 'IN STRATAER : year_emit_sai start/end', year_emit_sai_start, year_emit_sai_end
     189          WRITE(lunout, *) 'IN STRATAER : mth_emit_sai start/end', mth_emit_sai_start, mth_emit_sai_end
     190          WRITE(lunout, *) 'IN STRATAER : day_emit_sai start/end', day_emit_sai_start, day_emit_sai_end
     191          WRITE(lunout, *) 'IN STRATAER : xlat_min_sai', xlat_min_sai
     192          WRITE(lunout, *) 'IN STRATAER : xlat_max_sai', xlat_max_sai
     193          WRITE(lunout, *) 'IN STRATAER : xlon_sai', xlon_sai
     194        ENDIF
     195        IF(flag_emit == 4) THEN
     196          WRITE(lunout, *) 'IN STRATAER : nSpeciesErupt: ', nSpeciesErupt
     197          WRITE(lunout, *) 'IN STRATAER : id_species = ', id_species
     198          WRITE(lunout, *) 'IN STRATAER : id_HCl = ', id_HCl
     199          WRITE(lunout, *) 'IN STRATAER : id_HBr = ', id_HBr
     200          WRITE(lunout, *) 'IN STRATAER : id_NOx = ', id_NOx
     201          WRITE(lunout, *) 'IN STRATAER : id_H2O = ', id_H2O
     202          WRITE(lunout, *) 'IN STRATAER : m_Chlore_emiss_vol = ', m_Chlore_emiss_vol
     203          WRITE(lunout, *) 'IN STRATAER : m_Brome_emiss_vol = ', m_Brome_emiss_vol
     204          WRITE(lunout, *) 'IN STRATAER : m_NOx_emiss_vol = ', m_NOx_emiss_vol
     205          WRITE(lunout, *) 'IN STRATAER : m_H2O_emiss_vol = ', m_H2O_emiss_vol
     206        ENDIF
     207      endif
    207208    ENDIF ! if master
    208209
    209     WRITE(lunout,*) 'IN STRATAER_EMISS END'
     210    WRITE(lunout, *) 'IN STRATAER_EMISS END'
    210211  END SUBROUTINE strataer_emiss_init
    211212
    … …  
    222223
    223224    ! local var
    224     REAL                :: lat_reg_deg,lon_reg_deg ! lat and lon of grid points in degree
    225     INTEGER             :: ieru, i, j
     225    REAL :: lat_reg_deg, lon_reg_deg ! lat and lon of grid points in degree
     226    INTEGER :: ieru, i, j
    226227
    227228    ALLOCATE(ponde_lonlat_vol(nErupt))
    228229
    229230    !Compute lon/lat ponderation for injection
    230     dlat=180./2./FLOAT(nbp_lat)   ! d latitude in degree
    231     dlon=360./2./FLOAT(nbp_lon)   ! d longitude in degree
    232     WRITE(lunout,*) 'IN STRATAER_INIT dlat=',dlat,'dlon=',dlon
    233     WRITE(lunout,*) 'IN STRATAER_INIT nErupt=',nErupt
    234     WRITE(lunout,*) 'IN STRATAER_INIT xlat_min=',xlat_min_vol,'xlat_max=',xlat_max_vol
    235     WRITE(lunout,*) 'IN STRATAER_INIT xlon_min=',xlon_min_vol,'xlon_max=',xlon_max_vol
    236 
    237     DO ieru=1, nErupt
    238        ponde_lonlat_vol(ieru) = 0
    239        DO i=1,nbp_lon
    240           lon_reg_deg = lon_reg(i)*180./RPI
    241           DO j=1,nbp_lat
    242              lat_reg_deg = lat_reg(j)*180./RPI
    243              IF  ( lat_reg_deg>=xlat_min_vol(ieru)-dlat .AND. lat_reg_deg<xlat_max_vol(ieru)+dlat .AND. &
    244                   lon_reg_deg>=xlon_min_vol(ieru)-dlon .AND. lon_reg_deg<xlon_max_vol(ieru)+dlon ) THEN
    245                 ponde_lonlat_vol(ieru) = ponde_lonlat_vol(ieru) + 1
    246              ENDIF
    247           ENDDO
    248        ENDDO
    249        IF(ponde_lonlat_vol(ieru) == 0) THEN
    250           WRITE(lunout,*) 'STRATAER_INIT ERROR: no grid point found for eruption ieru=',ieru
    251        ENDIF
     231    dlat = 180. / 2. / FLOAT(nbp_lat)   ! d latitude in degree
     232    dlon = 360. / 2. / FLOAT(nbp_lon)   ! d longitude in degree
     233    WRITE(lunout, *) 'IN STRATAER_INIT dlat=', dlat, 'dlon=', dlon
     234    WRITE(lunout, *) 'IN STRATAER_INIT nErupt=', nErupt
     235    WRITE(lunout, *) 'IN STRATAER_INIT xlat_min=', xlat_min_vol, 'xlat_max=', xlat_max_vol
     236    WRITE(lunout, *) 'IN STRATAER_INIT xlon_min=', xlon_min_vol, 'xlon_max=', xlon_max_vol
     237
     238    DO ieru = 1, nErupt
     239      ponde_lonlat_vol(ieru) = 0
     240      DO i = 1, nbp_lon
     241        lon_reg_deg = lon_reg(i) * 180. / RPI
     242        DO j = 1, nbp_lat
     243          lat_reg_deg = lat_reg(j) * 180. / RPI
     244          IF  (lat_reg_deg>=xlat_min_vol(ieru) - dlat .AND. lat_reg_deg<xlat_max_vol(ieru) + dlat .AND. &
     245                  lon_reg_deg>=xlon_min_vol(ieru) - dlon .AND. lon_reg_deg<xlon_max_vol(ieru) + dlon) THEN
     246            ponde_lonlat_vol(ieru) = ponde_lonlat_vol(ieru) + 1
     247          ENDIF
     248        ENDDO
     249      ENDDO
     250      IF(ponde_lonlat_vol(ieru) == 0) THEN
     251        WRITE(lunout, *) 'STRATAER_INIT ERROR: no grid point found for eruption ieru=', ieru
     252      ENDIF
    252253    ENDDO !ieru
    253254
    254     WRITE(lunout,*) 'IN STRATAER_PONDE_INIT ponde_lonlat: ', ponde_lonlat_vol
     255    WRITE(lunout, *) 'IN STRATAER_PONDE_INIT ponde_lonlat: ', ponde_lonlat_vol
    255256
    256257  END SUBROUTINE strataer_ponde_init

  • LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/strataer_local_var_mod.F90

    r5110 r5111  
    11MODULE strataer_local_var_mod
    22  ! This module contains strato microphysic model parameters & variables
    3 
     3  USE lmdz_abort_physic, ONLY: abort_physic
    44  IMPLICIT NONE
    55

  • LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/traccoag_mod.F90

    r5110 r5111  
    11MODULE traccoag_mod
    22
    3 ! This module calculates the concentration of aerosol particles in certain size bins
    4 ! considering coagulation and sedimentation.
    5 
     3  ! This module calculates the concentration of aerosol particles in certain size bins
     4  ! considering coagulation and sedimentation.
     5  USE lmdz_abort_physic, ONLY: abort_physic
    66CONTAINS
    77
    88  SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, &
    9        presnivs, xlat, xlon, pphis, pphi, &
    10        t_seri, pplay, paprs, sh, rh, tr_seri)
    11    
    12     USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, & 
    13    budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part, &
    14    R2SO4B, DENSO4B, f_r_wetB, sulfmmr, SAD_sulfate, sulfmmr_mode, nd_mode
    15    
     9          presnivs, xlat, xlon, pphis, pphi, &
     10          t_seri, pplay, paprs, sh, rh, tr_seri)
     11
     12    USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, &
     13            budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part, &
     14            R2SO4B, DENSO4B, f_r_wetB, sulfmmr, SAD_sulfate, sulfmmr_mode, nd_mode
     15
    1616    USE dimphy
    1717    USE infotrac_phy, ONLY: nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat
    … …  
    1919    USE geometry_mod, ONLY: cell_area, boundslat
    2020    USE lmdz_grid_phy
    21     USE lmdz_phys_mpi_data, ONLY:  is_mpi_root
     21    USE lmdz_phys_mpi_data, ONLY: is_mpi_root
    2222    USE lmdz_phys_para, only: gather, scatter
    2323    USE phys_cal_mod, ONLY: year_len, year_cur, mth_cur, day_cur, hour
    … …  
    2727    USE print_control_mod, ONLY: lunout
    2828    USE strataer_local_var_mod
    29    
     29
    3030    IMPLICIT NONE
    3131
    32 ! Input argument
    33 !---------------
    34     REAL,INTENT(IN)    :: pdtphys    ! Pas d'integration pour la physique (seconde)
    35     REAL,INTENT(IN)    :: gmtime     ! Heure courante
    36     LOGICAL,INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
    37     INTEGER,INTENT(IN) :: julien     ! Jour julien
    38 
    39     REAL,DIMENSION(klev),INTENT(IN)        :: presnivs! pressions approximat. des milieux couches (en PA)
    40     REAL,DIMENSION(klon),INTENT(IN)        :: xlat    ! latitudes pour chaque point
    41     REAL,DIMENSION(klon),INTENT(IN)        :: xlon    ! longitudes pour chaque point
    42     REAL,DIMENSION(klon),INTENT(IN)        :: pphis   ! geopotentiel du sol
    43     REAL,DIMENSION(klon,klev),INTENT(IN)  :: pphi    ! geopotentiel de chaque couche
    44 
    45     REAL,DIMENSION(klon,klev),INTENT(IN)  :: t_seri  ! Temperature
    46     REAL,DIMENSION(klon,klev),INTENT(IN)  :: pplay   ! pression pour le mileu de chaque couche (en Pa)
    47     REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
    48     REAL,DIMENSION(klon,klev),INTENT(IN)  :: sh      ! humidite specifique
    49     REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative   
    50 
    51 ! Output argument
    52 !----------------
    53     REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA] 
    54 
    55 ! Local variables
    56 !----------------
    57     REAL                                   :: m_aer_emiss_vol_daily ! daily injection mass emission
    58     REAL                                   :: m_aer               ! aerosol mass
    59     INTEGER                                :: it, k, i, j, ilon, ilev, itime, i_int, ieru
    60     LOGICAL,DIMENSION(klon,klev)          :: is_strato           ! true = above tropopause, false = below
    61     REAL,DIMENSION(klon,klev)              :: m_air_gridbox       ! mass of air in every grid box [kg]
    62     REAL,DIMENSION(klon_glo,klev,nbtr)     :: tr_seri_glo         ! Concentration Traceur [U/KgA] 
    63     REAL,DIMENSION(klev+1)                :: altLMDz             ! altitude of layer interfaces in m
    64     REAL,DIMENSION(klev)                  :: f_lay_emiss         ! fraction of emission for every vertical layer
    65     REAL                                   :: f_lay_sum           ! sum of layer emission fractions
    66     REAL                                   :: alt                 ! altitude for integral calculation
    67     INTEGER,PARAMETER                      :: n_int_alt=10        ! number of subintervals for integration over Gaussian emission profile
    68     REAL,DIMENSION(nbtr_bin)              :: r_bin               ! particle radius in size bin [m]
    69     REAL,DIMENSION(nbtr_bin)              :: r_lower             ! particle radius at lower bin boundary [m]
    70     REAL,DIMENSION(nbtr_bin)              :: r_upper             ! particle radius at upper bin boundary [m]
    71     REAL,DIMENSION(nbtr_bin)              :: m_part_dry          ! mass of one dry particle in size bin [kg]
    72     REAL                                   :: zrho                ! Density of air [kg/m3]
    73     REAL                                   :: zdz                 ! thickness of atm. model layer in m
    74     REAL,DIMENSION(klev)                  :: zdm                 ! mass of atm. model layer in kg
    75     REAL,DIMENSION(klon,klev)              :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
    76     REAL                                   :: emission            ! emission
    77     REAL                                   :: theta_min, theta_max ! for SAI computation between two latitudes
    78     REAL                                   :: dlat_loc
    79     REAL                                   :: latmin,latmax,lonmin,lonmax ! lat/lon min/max for injection
    80     REAL                                   :: sigma_alt, altemiss ! injection altitude + sigma for distrib
    81     REAL                                   :: pdt,stretchlong     ! physic timestep, stretch emission over one day
    82    
    83     INTEGER                                :: injdur_sai          ! injection duration for SAI case [days]
    84     INTEGER                                :: yr, is_bissext
     32    ! Input argument
     33    !---------------
     34    REAL, INTENT(IN) :: pdtphys    ! Pas d'integration pour la physique (seconde)
     35    REAL, INTENT(IN) :: gmtime     ! Heure courante
     36    LOGICAL, INTENT(IN) :: debutphy   ! le flag de l'initialisation de la physique
     37    INTEGER, INTENT(IN) :: julien     ! Jour julien
     38
     39    REAL, DIMENSION(klev), INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA)
     40    REAL, DIMENSION(klon), INTENT(IN) :: xlat    ! latitudes pour chaque point
     41    REAL, DIMENSION(klon), INTENT(IN) :: xlon    ! longitudes pour chaque point
     42    REAL, DIMENSION(klon), INTENT(IN) :: pphis   ! geopotentiel du sol
     43    REAL, DIMENSION(klon, klev), INTENT(IN) :: pphi    ! geopotentiel de chaque couche
     44
     45    REAL, DIMENSION(klon, klev), INTENT(IN) :: t_seri  ! Temperature
     46    REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay   ! pression pour le mileu de chaque couche (en Pa)
     47    REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs   ! pression pour chaque inter-couche (en Pa)
     48    REAL, DIMENSION(klon, klev), INTENT(IN) :: sh      ! humidite specifique
     49    REAL, DIMENSION(klon, klev), INTENT(IN) :: rh      ! humidite relative
     50
     51    ! Output argument
     52    !----------------
     53    REAL, DIMENSION(klon, klev, nbtr), INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
     54
     55    ! Local variables
     56    !----------------
     57    REAL :: m_aer_emiss_vol_daily ! daily injection mass emission
     58    REAL :: m_aer               ! aerosol mass
     59    INTEGER :: it, k, i, j, ilon, ilev, itime, i_int, ieru
     60    LOGICAL, DIMENSION(klon, klev) :: is_strato           ! true = above tropopause, false = below
     61    REAL, DIMENSION(klon, klev) :: m_air_gridbox       ! mass of air in every grid box [kg]
     62    REAL, DIMENSION(klon_glo, klev, nbtr) :: tr_seri_glo         ! Concentration Traceur [U/KgA]
     63    REAL, DIMENSION(klev + 1) :: altLMDz             ! altitude of layer interfaces in m
     64    REAL, DIMENSION(klev) :: f_lay_emiss         ! fraction of emission for every vertical layer
     65    REAL :: f_lay_sum           ! sum of layer emission fractions
     66    REAL :: alt                 ! altitude for integral calculation
     67    INTEGER, PARAMETER :: n_int_alt = 10        ! number of subintervals for integration over Gaussian emission profile
     68    REAL, DIMENSION(nbtr_bin) :: r_bin               ! particle radius in size bin [m]
     69    REAL, DIMENSION(nbtr_bin) :: r_lower             ! particle radius at lower bin boundary [m]
     70    REAL, DIMENSION(nbtr_bin) :: r_upper             ! particle radius at upper bin boundary [m]
     71    REAL, DIMENSION(nbtr_bin) :: m_part_dry          ! mass of one dry particle in size bin [kg]
     72    REAL :: zrho                ! Density of air [kg/m3]
     73    REAL :: zdz                 ! thickness of atm. model layer in m
     74    REAL, DIMENSION(klev) :: zdm                 ! mass of atm. model layer in kg
     75    REAL, DIMENSION(klon, klev) :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
     76    REAL :: emission            ! emission
     77    REAL :: theta_min, theta_max ! for SAI computation between two latitudes
     78    REAL :: dlat_loc
     79    REAL :: latmin, latmax, lonmin, lonmax ! lat/lon min/max for injection
     80    REAL :: sigma_alt, altemiss ! injection altitude + sigma for distrib
     81    REAL :: pdt, stretchlong     ! physic timestep, stretch emission over one day
     82
     83    INTEGER :: injdur_sai          ! injection duration for SAI case [days]
     84    INTEGER :: yr, is_bissext
    8585
    8686    IF (is_mpi_root .AND. flag_verbose_strataer) THEN
    87        WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
    88        WRITE(lunout,*) 'IN traccoag flag_emit: ',flag_emit
     87      WRITE(lunout, *) 'in traccoag: date from phys_cal_mod =', year_cur, '-', mth_cur, '-', day_cur, '-', hour
     88      WRITE(lunout, *) 'IN traccoag flag_emit: ', flag_emit
    8989    ENDIF
    90    
     90
    9191    !   radius [m]
    92     DO it=1, nbtr_bin
    93       r_bin(it)=mdw(it)/2.
    94     ENDDO
    95 
    96 !--set boundaries of size bins
    97     DO it=1, nbtr_bin
    98     IF (it==1) THEN
    99       r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
    100       r_lower(it)=r_bin(it)**2./r_upper(it)
    101     ELSEIF (it==nbtr_bin) THEN
    102       r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
    103       r_upper(it)=r_bin(it)**2./r_lower(it)
     92    DO it = 1, nbtr_bin
     93      r_bin(it) = mdw(it) / 2.
     94    ENDDO
     95
     96    !--set boundaries of size bins
     97    DO it = 1, nbtr_bin
     98      IF (it==1) THEN
     99        r_upper(it) = sqrt(r_bin(it + 1) * r_bin(it))
     100        r_lower(it) = r_bin(it)**2. / r_upper(it)
     101      ELSEIF (it==nbtr_bin) THEN
     102        r_lower(it) = sqrt(r_bin(it) * r_bin(it - 1))
     103        r_upper(it) = r_bin(it)**2. / r_lower(it)
     104      ELSE
     105        r_lower(it) = sqrt(r_bin(it) * r_bin(it - 1))
     106        r_upper(it) = sqrt(r_bin(it + 1) * r_bin(it))
     107      ENDIF
     108    ENDDO
     109
     110    IF (debutphy .and. is_mpi_root) THEN
     111      DO it = 1, nbtr_bin
     112        WRITE(lunout, *) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'
     113      ENDDO
     114    ENDIF
     115
     116    !--initialising logical is_strato from stratomask
     117    is_strato(:, :) = .FALSE.
     118    WHERE (stratomask>0.5) is_strato = .TRUE.
     119
     120    IF(flag_new_strat_compo) THEN
     121      IF(debutphy) WRITE(lunout, *) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_new_strat_compo
     122      ! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc...
     123      CALL stracomp_kelvin(sh, t_seri, pplay)
    104124    ELSE
    105       r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
    106       r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
     125      IF(debutphy) WRITE(lunout, *) 'traccoag: COMPO from Bekki 2D model', flag_new_strat_compo
     126      ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
     127      ! H2SO4 mass fraction in aerosol (%)
     128      CALL stracomp(sh, t_seri, pplay)
     129
     130      ! aerosol density (gr/cm3)
     131      CALL denh2sa(t_seri)
     132
     133      ! compute factor for converting dry to wet radius (for every grid box)
     134      f_r_wet(:, :) = (dens_aer_dry / (DENSO4(:, :) * 1000.) / (R2SO4(:, :) / 100.))**(1. / 3.)
    107135    ENDIF
    108     ENDDO
    109 
    110     IF (debutphy .and. is_mpi_root) THEN
    111       DO it=1, nbtr_bin
    112         WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from',  r_lower(it), 'to', r_upper(it), ')'
    113       ENDDO
    114     ENDIF
    115 
    116 !--initialising logical is_strato from stratomask
    117     is_strato(:,:)=.FALSE.
    118     WHERE (stratomask>0.5) is_strato=.TRUE.
    119 
    120     IF(flag_new_strat_compo) THEN
    121        IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_new_strat_compo
    122        ! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc...
    123        CALL stracomp_kelvin(sh,t_seri,pplay)
    124     ELSE
    125        IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO from Bekki 2D model', flag_new_strat_compo
    126        ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
    127        ! H2SO4 mass fraction in aerosol (%)
    128        CALL stracomp(sh,t_seri,pplay)
    129        
    130        ! aerosol density (gr/cm3)
    131        CALL denh2sa(t_seri)
    132        
    133        ! compute factor for converting dry to wet radius (for every grid box)
    134        f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.)
    135     ENDIF
    136    
    137 !--calculate mass of air in every grid box
    138     DO ilon=1, klon
    139        DO ilev=1, klev
    140           m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
    141        ENDDO
    142     ENDDO
    143    
    144 !--initialise emission diagnostics
    145     if (nErupt > 0 .and. (flag_emit == 1 .or. flag_emit == 4)) budg_emi(:,1)=0.0
    146     budg_emi_ocs(:)=0.0
    147     budg_emi_so2(:)=0.0
    148     budg_emi_h2so4(:)=0.0
    149     budg_emi_part(:)=0.0
    150 
    151 !--sulfur emission, depending on chosen scenario (flag_emit)
     136
     137    !--calculate mass of air in every grid box
     138    DO ilon = 1, klon
     139      DO ilev = 1, klev
     140        m_air_gridbox(ilon, ilev) = (paprs(ilon, ilev) - paprs(ilon, ilev + 1)) / RG * cell_area(ilon)
     141      ENDDO
     142    ENDDO
     143
     144    !--initialise emission diagnostics
     145    if (nErupt > 0 .and. (flag_emit == 1 .or. flag_emit == 4)) budg_emi(:, 1) = 0.0
     146    budg_emi_ocs(:) = 0.0
     147    budg_emi_so2(:) = 0.0
     148    budg_emi_h2so4(:) = 0.0
     149    budg_emi_part(:) = 0.0
     150
     151    !--sulfur emission, depending on chosen scenario (flag_emit)
    152152    SELECT CASE(flag_emit)
    153153
    … …  
    158158      !--only emit on day of eruption
    159159      ! stretch emission over one day of Pinatubo eruption
    160        DO ieru=1, nErupt
    161           IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
    162                day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN
    163 
    164              ! daily injection mass emission
    165              m_aer=m_aer_emiss_vol(ieru,1)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
    166              !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
    167              m_aer=m_aer*(mSO2mol/mSatom)
    168              
    169              WRITE(lunout,*) 'IN traccoag m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru,1), &
    170                   'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', &
    171                   (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer,'ieru=',ieru
    172              WRITE(lunout,*) 'IN traccoag, dlon=',dlon
    173              
    174              latmin=xlat_min_vol(ieru)
    175              latmax=xlat_max_vol(ieru)
    176              lonmin=xlon_min_vol(ieru)
    177              lonmax=xlon_max_vol(ieru)
    178              altemiss = altemiss_vol(ieru)
    179              sigma_alt = sigma_alt_vol(ieru)
    180              pdt=pdtphys
    181              ! stretch emission over one day of eruption
    182              stretchlong = 1.
    183              
    184              CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,tr_seri,&
    185                   m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
    186                   stretchlong,1,0)
    187              
    188           ENDIF ! emission period
    189        ENDDO ! eruption number
    190        
     160      DO ieru = 1, nErupt
     161        IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
     162                day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru) + injdur)) THEN
     163
     164          ! daily injection mass emission
     165          m_aer = m_aer_emiss_vol(ieru, 1) / (REAL(injdur) * REAL(ponde_lonlat_vol(ieru)))
     166          !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
     167          m_aer = m_aer * (mSO2mol / mSatom)
     168
     169          WRITE(lunout, *) 'IN traccoag m_aer_emiss_vol(ieru)=', m_aer_emiss_vol(ieru, 1), &
     170                  'ponde_lonlat_vol(ieru)=', ponde_lonlat_vol(ieru), '(injdur*ponde_lonlat_vol(ieru))', &
     171                  (injdur * ponde_lonlat_vol(ieru)), 'm_aer_emiss_vol_daily=', m_aer, 'ieru=', ieru
     172          WRITE(lunout, *) 'IN traccoag, dlon=', dlon
     173
     174          latmin = xlat_min_vol(ieru)
     175          latmax = xlat_max_vol(ieru)
     176          lonmin = xlon_min_vol(ieru)
     177          lonmax = xlon_max_vol(ieru)
     178          altemiss = altemiss_vol(ieru)
     179          sigma_alt = sigma_alt_vol(ieru)
     180          pdt = pdtphys
     181          ! stretch emission over one day of eruption
     182          stretchlong = 1.
     183
     184          CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, tr_seri, &
     185                  m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, &
     186                  stretchlong, 1, 0)
     187
     188        ENDIF ! emission period
     189      ENDDO ! eruption number
     190
    191191    CASE(2) ! stratospheric aerosol injections (SAI)
    192192
    193      ! Computing duration of SAI in days...
    194      ! ... starting from 0...
    195      injdur_sai = 0
    196      ! ... then adding whole years from first to (n-1)th...
    197      DO yr = year_emit_sai_start, year_emit_sai_end-1
    198        ! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0)
    199        is_bissext = (MOD(yr,4)==0) .AND. (MOD(yr,100) /= 0 .OR. MOD(yr,400) == 0)
    200        injdur_sai = injdur_sai+365+is_bissext
    201      ENDDO
    202      ! ... then subtracting part of the first year where no injection yet...
    203      is_bissext = (MOD(year_emit_sai_start,4)==0) .AND. (MOD(year_emit_sai_start,100) /= 0 .OR. MOD(year_emit_sai_start,400) == 0)
    204      SELECT CASE(mth_emit_sai_start)
    205      CASE(2)
    206         injdur_sai = injdur_sai-31
    207      CASE(3)
    208         injdur_sai = injdur_sai-31-28-is_bissext
    209      CASE(4)
    210         injdur_sai = injdur_sai-31-28-is_bissext-31
    211      CASE(5)
    212         injdur_sai = injdur_sai-31-28-is_bissext-31-30
    213      CASE(6)
    214         injdur_sai = injdur_sai-31-28-is_bissext-31-30-31
    215      CASE(7)
    216         injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30
    217      CASE(8)
    218         injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31
    219      CASE(9)
    220         injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31
    221      CASE(10)
    222         injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30
    223      CASE(11)
    224         injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31
    225      CASE(12)
    226         injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31-30
    227      END SELECT
    228      injdur_sai = injdur_sai-day_emit_sai_start+1
    229      ! ... then adding part of the n-th year
    230      is_bissext = (MOD(year_emit_sai_end,4)==0) .AND. (MOD(year_emit_sai_end,100) /= 0 .OR. MOD(year_emit_sai_end,400) == 0)
    231      SELECT CASE(mth_emit_sai_end)
    232      CASE(2)
    233         injdur_sai = injdur_sai+31
    234      CASE(3)
    235         injdur_sai = injdur_sai+31+28+is_bissext
    236      CASE(4)
    237         injdur_sai = injdur_sai+31+28+is_bissext+31
    238      CASE(5)
    239         injdur_sai = injdur_sai+31+28+is_bissext+31+30
    240      CASE(6)
    241         injdur_sai = injdur_sai+31+28+is_bissext+31+30+31
    242      CASE(7)
    243         injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30
    244      CASE(8)
    245         injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31
    246      CASE(9)
    247         injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31
    248      CASE(10)
    249         injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30
    250      CASE(11)
    251         injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31
    252      CASE(12)
    253         injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31+30
    254      END SELECT
    255      injdur_sai = injdur_sai+day_emit_sai_end
    256      ! A security: are SAI dates of injection consistent?
    257      IF (injdur_sai <= 0) THEN
    258         CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.',1)
    259      ENDIF
    260      ! Injection in itself
    261      IF (( year_emit_sai_start <= year_cur ) &
    262         .AND. ( year_cur <= year_emit_sai_end ) &
    263         .AND. ( mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur ) &
    264         .AND. ( mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end ) &
    265         .AND. ( day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur ) &
    266         .AND. ( day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end )) THEN
    267        
    268        m_aer=m_aer_emiss_sai
    269        !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
    270        m_aer=m_aer*(mSO2mol/mSatom)
    271        
    272        latmin=xlat_sai
    273        latmax=xlat_sai
    274        lonmin=xlon_sai
    275        lonmax=xlon_sai
    276        altemiss = altemiss_sai
    277        sigma_alt = sigma_alt_sai
    278        pdt=0.
    279        ! stretch emission over whole year (360d)
    280        stretchlong=FLOAT(year_len)
    281        
    282        CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
    283             m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
    284             stretchlong,1,0)
    285        
    286        budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
    287      ENDIF ! Condition over injection dates
     193      ! Computing duration of SAI in days...
     194      ! ... starting from 0...
     195      injdur_sai = 0
     196      ! ... then adding whole years from first to (n-1)th...
     197      DO yr = year_emit_sai_start, year_emit_sai_end - 1
     198        ! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0)
     199        is_bissext = (MOD(yr, 4)==0) .AND. (MOD(yr, 100) /= 0 .OR. MOD(yr, 400) == 0)
     200        injdur_sai = injdur_sai + 365 + is_bissext
     201      ENDDO
     202      ! ... then subtracting part of the first year where no injection yet...
     203      is_bissext = (MOD(year_emit_sai_start, 4)==0) .AND. (MOD(year_emit_sai_start, 100) /= 0 .OR. MOD(year_emit_sai_start, 400) == 0)
     204      SELECT CASE(mth_emit_sai_start)
     205      CASE(2)
     206        injdur_sai = injdur_sai - 31
     207      CASE(3)
     208        injdur_sai = injdur_sai - 31 - 28 - is_bissext
     209      CASE(4)
     210        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31
     211      CASE(5)
     212        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30
     213      CASE(6)
     214        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31
     215      CASE(7)
     216        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30
     217      CASE(8)
     218        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31
     219      CASE(9)
     220        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31
     221      CASE(10)
     222        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30
     223      CASE(11)
     224        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 - 31
     225      CASE(12)
     226        injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 - 31 - 30
     227      END SELECT
     228      injdur_sai = injdur_sai - day_emit_sai_start + 1
     229      ! ... then adding part of the n-th year
     230      is_bissext = (MOD(year_emit_sai_end, 4)==0) .AND. (MOD(year_emit_sai_end, 100) /= 0 .OR. MOD(year_emit_sai_end, 400) == 0)
     231      SELECT CASE(mth_emit_sai_end)
     232      CASE(2)
     233        injdur_sai = injdur_sai + 31
     234      CASE(3)
     235        injdur_sai = injdur_sai + 31 + 28 + is_bissext
     236      CASE(4)
     237        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31
     238      CASE(5)
     239        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30
     240      CASE(6)
     241        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31
     242      CASE(7)
     243        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30
     244      CASE(8)
     245        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31
     246      CASE(9)
     247        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31
     248      CASE(10)
     249        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30
     250      CASE(11)
     251        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31
     252      CASE(12)
     253        injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30
     254      END SELECT
     255      injdur_sai = injdur_sai + day_emit_sai_end
     256      ! A security: are SAI dates of injection consistent?
     257      IF (injdur_sai <= 0) THEN
     258        CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.', 1)
     259      ENDIF
     260      ! Injection in itself
     261      IF ((year_emit_sai_start <= year_cur) &
     262              .AND. (year_cur <= year_emit_sai_end) &
     263              .AND. (mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur) &
     264              .AND. (mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end) &
     265              .AND. (day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur) &
     266              .AND. (day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end)) THEN
     267
     268        m_aer = m_aer_emiss_sai
     269        !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
     270        m_aer = m_aer * (mSO2mol / mSatom)
     271
     272        latmin = xlat_sai
     273        latmax = xlat_sai
     274        lonmin = xlon_sai
     275        lonmax = xlon_sai
     276        altemiss = altemiss_sai
     277        sigma_alt = sigma_alt_sai
     278        pdt = 0.
     279        ! stretch emission over whole year (360d)
     280        stretchlong = FLOAT(year_len)
     281
     282        CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, m_air_gridbox, tr_seri, &
     283                m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, &
     284                stretchlong, 1, 0)
     285
     286        budg_emi_so2(:) = budg_emi(:, 1) * mSatom / mSO2mol
     287      ENDIF ! Condition over injection dates
    288288
    289289    CASE(3) ! --- SAI injection over a single band of longitude and between
    290             !     lat_min and lat_max
    291 
    292        m_aer=m_aer_emiss_sai
    293        !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
    294        m_aer=m_aer*(mSO2mol/mSatom)
    295 
    296        latmin=xlat_min_sai
    297        latmax=xlat_max_sai
    298        lonmin=xlon_sai
    299        lonmax=xlon_sai
    300        altemiss = altemiss_sai
    301        sigma_alt = sigma_alt_sai
    302        pdt=0.
    303        ! stretch emission over whole year (360d)
    304        stretchlong=FLOAT(year_len)
    305 
    306        CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
    307             m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
    308             stretchlong,1,0)
    309 
    310        budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
    311        
     290      !     lat_min and lat_max
     291
     292      m_aer = m_aer_emiss_sai
     293      !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
     294      m_aer = m_aer * (mSO2mol / mSatom)
     295
     296      latmin = xlat_min_sai
     297      latmax = xlat_max_sai
     298      lonmin = xlon_sai
     299      lonmax = xlon_sai
     300      altemiss = altemiss_sai
     301      sigma_alt = sigma_alt_sai
     302      pdt = 0.
     303      ! stretch emission over whole year (360d)
     304      stretchlong = FLOAT(year_len)
     305
     306      CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, m_air_gridbox, tr_seri, &
     307              m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, &
     308              stretchlong, 1, 0)
     309
     310      budg_emi_so2(:) = budg_emi(:, 1) * mSatom / mSO2mol
     311
    312312    END SELECT ! emission scenario (flag_emit)
    313313
    314 !--read background concentrations of OCS and SO2 and lifetimes from input file
    315 !--update the variables defined in phys_local_var_mod
    316     CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
    317 
    318 !--convert OCS to SO2 in the stratosphere
    319     CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
    320 
    321 !--convert SO2 to H2SO4
    322     CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
    323 
    324 !--common routine for nucleation and condensation/evaporation with adaptive timestep
    325     CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
    326 
    327 !--CALL coagulation routine
    328     CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
    329 
    330 !--CALL sedimentation routine
     314    !--read background concentrations of OCS and SO2 and lifetimes from input file
     315    !--update the variables defined in phys_local_var_mod
     316    CALL interp_sulf_input(debutphy, pdtphys, paprs, tr_seri)
     317
     318    !--convert OCS to SO2 in the stratosphere
     319    CALL ocs_to_so2(pdtphys, tr_seri, t_seri, pplay, paprs, is_strato)
     320
     321    !--convert SO2 to H2SO4
     322    CALL so2_to_h2so4(pdtphys, tr_seri, t_seri, pplay, paprs, is_strato)
     323
     324    !--common routine for nucleation and condensation/evaporation with adaptive timestep
     325    CALL micphy_tstep(pdtphys, tr_seri, t_seri, pplay, paprs, rh, is_strato)
     326
     327    !--CALL coagulation routine
     328    CALL coagulate(pdtphys, mdw, tr_seri, t_seri, pplay, dens_aer, is_strato)
     329
     330    !--CALL sedimentation routine
    331331    CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
    332332
    333 !--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
    334     surf_PM25_sulf(:)=0.0
    335     DO i=1,klon
    336       DO it=1, nbtr_bin
     333    !--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
     334    surf_PM25_sulf(:) = 0.0
     335    DO i = 1, klon
     336      DO it = 1, nbtr_bin
    337337        IF (mdw(it) < 2.5e-6) THEN
    338338          !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
    339339          !assume that particles consist of ammonium sulfate at the surface (132g/mol)
    340340          !and are dry at T = 20 deg. C and 50 perc. humidity
    341           surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
    342    *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 &
    343    *pplay(i,1)/t_seri(i,1)/RD*1.e9
     341          surf_PM25_sulf(i) = surf_PM25_sulf(i) + tr_seri(i, 1, it + nbtr_sulgas) &
     342                  * 132. / 98. * dens_aer_dry * 4. / 3. * RPI * (mdw(it) / 2.)**3 &
     343                  * pplay(i, 1) / t_seri(i, 1) / RD * 1.e9
    344344        ENDIF
    345345      ENDDO
    346346    ENDDO
    347    
    348 !--compute
    349 !     sulfmmr: Sulfate aerosol concentration (dry mixing ratio) (condensed H2SO4 mmr)
    350 !     SAD_sulfate: SAD all aerosols (cm2/cm3) (must be WET)
    351 !     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (ambiguous but based on sulfmmr, it mus be DRY(?) mmr)
    352 !     nd_mode: DRY(?) particle concentration in different modes (part/m3)
    353      sulfmmr(:,:)=0.0
    354      SAD_sulfate(:,:)=0.0
    355      sulfmmr_mode(:,:,:)=0.0
    356      nd_mode(:,:,:)=0.0
    357      
    358      DO i=1,klon
    359         DO j=1,klev
    360            DO it=1, nbtr_bin
    361               !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
    362               !assume that particles consist of ammonium sulfate at the surface (132g/mol)
    363               !and are dry at T = 20 deg. C and 50 perc. humidity
    364              
    365               !     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (based on sulfmmr, it must be DRY mmr)
    366               !     equivalent to condensed H2SO4 mmr= H2SO4 kg / kgA in bin it
    367               sulfmmr_mode(i,j,it) = tr_seri(i,j,it+nbtr_sulgas) &        ! [DRY part/kgA in bin it]
    368     *(4./3.)*RPI*(mdw(it)/2.)**3.   &                   ! [mdw: dry diameter in m]
    369     *dens_aer_dry                                       ! [dry aerosol mass density in kg/m3]
    370              
    371               !     sulfmmr: Sulfate aerosol concentration (dry mass mixing ratio)
    372               !     equivalent to total condensed H2SO4 mmr (H2SO4 kg / kgA
    373               sulfmmr(i,j) = sulfmmr(i,j) + sulfmmr_mode(i,j,it)
    374              
    375               !     nd_mode: particle concentration in different modes (DRY part/m3)
    376               nd_mode(i,j,it) = tr_seri(i,j,it+nbtr_sulgas) &             ! [DRY part/kgA in bin it]
    377    *pplay(i,j)/t_seri(i,j)/RD                           ! [air mass concentration in kg air /m3A]
    378              
    379               IF(flag_new_strat_compo) THEN
    380                  !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
    381                  SAD_sulfate(i,j) = SAD_sulfate(i,j) + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
    382     *4.*RPI*( mdw(it)*f_r_wetB(i,j,it)/2. )**2. &       ! [WET SA of part it in m2]
    383     *1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
    384               ELSE
    385                  !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
    386                  SAD_sulfate(i,j) = SAD_sulfate(i,j) + nd_mode(i,j,it) &     ! [DRY part/m3A (in bin it)]
    387     *4.*RPI*( mdw(it)*f_r_wet(i,j)/2. )**2. &           ! [WET SA of part it in m2]
    388     *1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
    389               ENDIF
    390            ENDDO
     347
     348    !--compute
     349    !     sulfmmr: Sulfate aerosol concentration (dry mixing ratio) (condensed H2SO4 mmr)
     350    !     SAD_sulfate: SAD all aerosols (cm2/cm3) (must be WET)
     351    !     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (ambiguous but based on sulfmmr, it mus be DRY(?) mmr)
     352    !     nd_mode: DRY(?) particle concentration in different modes (part/m3)
     353    sulfmmr(:, :) = 0.0
     354    SAD_sulfate(:, :) = 0.0
     355    sulfmmr_mode(:, :, :) = 0.0
     356    nd_mode(:, :, :) = 0.0
     357
     358    DO i = 1, klon
     359      DO j = 1, klev
     360        DO it = 1, nbtr_bin
     361          !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
     362          !assume that particles consist of ammonium sulfate at the surface (132g/mol)
     363          !and are dry at T = 20 deg. C and 50 perc. humidity
     364
     365          !     sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (based on sulfmmr, it must be DRY mmr)
     366          !     equivalent to condensed H2SO4 mmr= H2SO4 kg / kgA in bin it
     367          sulfmmr_mode(i, j, it) = tr_seri(i, j, it + nbtr_sulgas) &        ! [DRY part/kgA in bin it]
     368                  * (4. / 3.) * RPI * (mdw(it) / 2.)**3.   &                   ! [mdw: dry diameter in m]
     369                  * dens_aer_dry                                       ! [dry aerosol mass density in kg/m3]
     370
     371          !     sulfmmr: Sulfate aerosol concentration (dry mass mixing ratio)
     372          !     equivalent to total condensed H2SO4 mmr (H2SO4 kg / kgA
     373          sulfmmr(i, j) = sulfmmr(i, j) + sulfmmr_mode(i, j, it)
     374
     375          !     nd_mode: particle concentration in different modes (DRY part/m3)
     376          nd_mode(i, j, it) = tr_seri(i, j, it + nbtr_sulgas) &             ! [DRY part/kgA in bin it]
     377                  * pplay(i, j) / t_seri(i, j) / RD                           ! [air mass concentration in kg air /m3A]
     378
     379          IF(flag_new_strat_compo) THEN
     380            !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
     381            SAD_sulfate(i, j) = SAD_sulfate(i, j) + nd_mode(i, j, it) &     ! [DRY part/m3A (in bin it)]
     382                    * 4. * RPI * (mdw(it) * f_r_wetB(i, j, it) / 2.)**2. &       ! [WET SA of part it in m2]
     383                    * 1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
     384          ELSE
     385            !     SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
     386            SAD_sulfate(i, j) = SAD_sulfate(i, j) + nd_mode(i, j, it) &     ! [DRY part/m3A (in bin it)]
     387                    * 4. * RPI * (mdw(it) * f_r_wet(i, j) / 2.)**2. &           ! [WET SA of part it in m2]
     388                    * 1.e-2                                              ! conversion from m2/m3 to cm2/cm3A
     389          ENDIF
    391390        ENDDO
    392      ENDDO
    393      
     391      ENDDO
     392    ENDDO
     393
    394394  END SUBROUTINE traccoag
    395395
Note: See TracChangeset for help on using the changeset viewer.