Ignore:
Timestamp:
Oct 30, 2023, 5:37:00 PM (8 months ago)
Author:
Laurent Fairhead
Message:

Merge of ACC branch with 4740 revision from trunk

Location:
LMDZ6/branches/Portage_acc
Files:
1 deleted
8 edited
6 copied

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Portage_acc

  • LMDZ6/branches/Portage_acc/libf/phylmd/StratAer/aerophys.F90

    r3526 r4743  
    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/branches/Portage_acc/libf/phylmd/StratAer/interp_sulf_input.F90

    r3677 r4743  
    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
    … …  
    4243  REAL paprs_glo(klon_glo,klev+1)
    4344
    44   REAL, POINTER:: latitude(:)
     45  REAL, ALLOCATABLE :: latitude(:)
    4546! (of input data sorted in strictly ascending order)
    4647
    47   REAL, POINTER:: longitude(:)
     48  REAL, ALLOCATABLE :: longitude(:)
    4849! (of input data sorted in strictly ascending order)
    4950
    50   REAL, POINTER:: time(:)
     51  REAL, ALLOCATABLE :: time(:)
    5152! (of input data sorted in strictly ascending order)
    5253
    53   REAL, POINTER:: lev(:)
     54  REAL, ALLOCATABLE :: lev(:)
    5455! levels of input data
    5556
    … …  
    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/branches/Portage_acc/libf/phylmd/StratAer/micphy_tstep.F90

    r4446 r4743  
    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/branches/Portage_acc/libf/phylmd/StratAer/nucleation_tstep_mod.F90

    r4446 r4743  
    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/branches/Portage_acc/libf/phylmd/StratAer/ocs_to_so2.F90

    r3677 r4743  
    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/branches/Portage_acc/libf/phylmd/StratAer/so2_to_h2so4.F90

    r3677 r4743  
    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/branches/Portage_acc/libf/phylmd/StratAer/traccoag_mod.F90

    r4584 r4743  
    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
Note: See TracChangeset for help on using the changeset viewer.