Ignore:
Timestamp:
Nov 20, 2016, 2:15:32 PM (8 years ago)
Author:
oboucher
Message:

This revision concerns the StratAer? module and should not impact the rest of LMDz
Bug correction in interp_sulf_input.F90
Update of miecalc_aer.F90 and traccoag_mod.F90
Phytrac tracers are now dealt with in XIOS through the Fortran interface with minimal input in the xml
Making tracer groups in DefLists? for StratAer?

Location:
LMDZ5/trunk/libf/phylmd/StratAer
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r2695 r2704  
    7171!$OMP THREADPRIVATE(OCS_clim,SO2_clim)
    7272!
    73   INTEGER i, k, kk, ilon, ilev, j
     73  INTEGER i, k, kk, j
    7474  REAL p_bound
    7575
    … …  
    8989  5.53084382e+01,   3.35462635e+01,   0.0           /)
    9090!
    91   REAL, PARAMETER :: epsilon_OCS=1.0e-20     ! minimum OCS concentration [kg/kgA] for weighting of lifetime
    92   REAL, PARAMETER :: epsilon_SO2=1.0e-20     ! minimum SO2 concentration [kg/kgA] for weighting of lifetime
    93   REAL, PARAMETER :: min_OCS_lifetime= 3600. !minimum OCS lifetime [sec]
    94   REAL, PARAMETER :: min_SO2_lifetime=86400. !minimum SO2 lifetime [sec]
    95 
    9691 IF (.NOT.ALLOCATED(OCS_clim)) ALLOCATE(OCS_clim(klon,klev))
    9792 IF (.NOT.ALLOCATED(SO2_clim)) ALLOCATE(SO2_clim(klon,klev))
    … …  
    183178
    184179  !---regrid weighted lifetime and climatologies
    185   DO i=1, klon
     180  DO i=1, klon_glo
    186181    DO k=1, klev
    187182     OCS_lifetime_glo(i,k)=0.0
    … …  
    206201  ENDDO
    207202
    208   ENDIF ! is_mpi_root
     203  ENDIF !--is_mpi_root and is_omp_root
    209204
    210205!--keep memory of previous month
    … …  
    255250
    256251  !convert SO2_backgr_tend from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic
    257   DO ilon=1, klon
    258     DO ilev=1, klev
    259       SO2_backgr_tend(ilon,ilev)=SO2_backgr_tend(ilon,ilev)*mSatom/mSO2mol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
    260       OCS_backgr_tend(ilon,ilev)=OCS_backgr_tend(ilon,ilev)*mSatom/mOCSmol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys
     252  DO i=1, klon
     253    DO k=1, klev
     254      SO2_backgr_tend(i,k)=SO2_backgr_tend(i,k)*mSatom/mSO2mol*(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
     255      OCS_backgr_tend(i,k)=OCS_backgr_tend(i,k)*mSatom/mOCSmol*(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
    261256    ENDDO
    262257  ENDDO

  • LMDZ5/trunk/libf/phylmd/StratAer/miecalc_aer.F90

    r2690 r2704  
    2424  USE mod_grid_phy_lmdz, ONLY : klon_glo
    2525  USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
     26  USE print_control_mod, ONLY: prt_level, lunout
    2627
    2728  IMPLICIT NONE
    … …  
    7071  PARAMETER (Nbin=10)
    7172
    72 !--has to be consistent with dataset
    73   INTEGER nb_lambda_h2so4
    74   PARAMETER (nb_lambda_h2so4=62) !61, 235
    7573
    7674!---wavelengths STREAMER
    … …  
    131129  REAL final_w(1:NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW)
    132130
    133   INTEGER band, bandSW, bandLW
     131  INTEGER band, bandSW, bandLW, wavenumber
    134132
    135133!---wavelengths SW RRTM
    … …  
    156154  REAL ss_g(nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw)
    157155
    158   REAL wavenumber
    159   REAL, DIMENSION(nb_lambda_h2so4,4) :: ref_ind
    160 
     156  INTEGER, PARAMETER :: nb_lambda_h2so4=62
     157  REAL, DIMENSION (nb_lambda_h2so4,4) :: ref_ind
     158  !-- fichier h2so4_0.75_300.00_hummel_1988_p_q.dat
     159  ! -- wavenumber (cm-1), wavelength (um), n_r, n_i
     160  DATA ref_ind /                                &
     161   200.000,   50.0000,   2.01000,   6.5000E-01, &
     162   250.000,   40.0000,   1.94000,   6.3000E-01, &
     163   285.714,   35.0000,   1.72000,   5.2000E-01, &
     164   333.333,   30.0000,   1.73000,   2.9000E-01, &
     165   358.423,   27.9000,   1.78000,   2.5000E-01, &
     166   400.000,   25.0000,   1.84000,   2.4000E-01, &
     167   444.444,   22.5000,   1.82000,   2.9000E-01, &
     168   469.484,   21.3000,   1.79000,   2.5000E-01, &
     169   500.000,   20.0000,   1.81000,   2.3000E-01, &
     170   540.541,   18.5000,   1.92700,   3.0200E-01, &
     171   555.556,   18.0000,   1.95000,   4.1000E-01, &
     172   581.395,   17.2000,   1.72400,   5.9000E-01, &
     173   609.756,   16.4000,   1.52000,   4.1400E-01, &
     174   666.667,   15.0000,   1.59000,   2.1100E-01, &
     175   675.676,   14.8000,   1.61000,   2.0500E-01, &
     176   714.286,   14.0000,   1.64000,   1.9500E-01, &
     177   769.231,   13.0000,   1.69000,   1.9500E-01, &
     178   800.000,   12.5000,   1.74000,   1.9800E-01, &
     179   869.565,   11.5000,   1.89000,   3.7400E-01, &
     180   909.091,   11.0000,   1.67000,   4.8500E-01, &
     181   944.198,   10.5910,   1.72000,   3.4000E-01, &
     182  1000.000,   10.0000,   1.89000,   4.5500E-01, &
     183  1020.408,    9.8000,   1.91000,   6.8000E-01, &
     184  1052.632,    9.5000,   1.67000,   7.5000E-01, &
     185  1086.957,    9.2000,   1.60000,   5.8600E-01, &
     186  1111.111,    9.0000,   1.65000,   6.3300E-01, &
     187  1149.425,    8.7000,   1.53000,   7.7200E-01, &
     188  1176.471,    8.5000,   1.37000,   7.5500E-01, &
     189  1219.512,    8.2000,   1.20000,   6.4500E-01, &
     190  1265.823,    7.9000,   1.14000,   4.8800E-01, &
     191  1388.889,    7.2000,   1.21000,   1.7600E-01, &
     192  1538.462,    6.5000,   1.37000,   1.2800E-01, &
     193  1612.903,    6.2000,   1.42400,   1.6500E-01, &
     194  1666.667,    6.0000,   1.42500,   1.9500E-01, &
     195  1818.182,    5.5000,   1.33700,   1.8300E-01, &
     196  2000.000,    5.0000,   1.36000,   1.2100E-01, &
     197  2222.222,    4.5000,   1.38500,   1.2000E-01, &
     198  2500.000,    4.0000,   1.39800,   1.2600E-01, &
     199  2666.667,    3.7500,   1.39600,   1.3100E-01, &
     200  2857.143,    3.5000,   1.37600,   1.5800E-01, &
     201  2948.113,    3.3920,   1.35200,   1.5900E-01, &
     202  3125.000,    3.2000,   1.31100,   1.3500E-01, &
     203  3333.333,    3.0000,   1.29300,   9.5500E-02, &
     204  3703.704,    2.7000,   1.30300,   5.7000E-03, &
     205  4000.000,    2.5000,   1.34400,   3.7600E-03, &
     206  4444.444,    2.2500,   1.37000,   1.8000E-03, &
     207  5000.000,    2.0000,   1.38400,   1.2600E-03, &
     208  5555.556,    1.8000,   1.39000,   5.5000E-04, &
     209  6510.417,    1.5360,   1.40300,   1.3700E-04, &
     210  7692.308,    1.3000,   1.41000,   1.0000E-05, &
     211  9433.962,    1.0600,   1.42000,   1.5000E-06, &
     212 11627.907,    0.8600,   1.42500,   1.7900E-07, &
     213 14409.222,    0.6940,   1.42800,   1.9900E-08, &
     214 15797.788,    0.6330,   1.42900,   1.4700E-08, &
     215 18181.818,    0.5500,   1.43000,   1.0000E-08, &
     216 19417.476,    0.5150,   1.43100,   1.0000E-08, &
     217 20491.803,    0.4880,   1.43200,   1.0000E-08, &
     218 25000.000,    0.4000,   1.44000,   1.0000E-08, &
     219 29673.591,    0.3370,   1.45900,   1.0000E-08, &
     220 33333.333,    0.3000,   1.46900,   1.0000E-08, &
     221 40000.000,    0.2500,   1.48400,   1.0000E-08, &
     222 50000.000,    0.2000,   1.49800,   1.0000E-08 /
    161223!---------------------------------------------------------
    162224
    163225  IF (debut) THEN   
     226
    164227  !--initialising dry diameters to geometrically spaced mass/volume (see Jacobson 1994)
    165228      mdw(1)=mdwmin
    … …  
    176239      PRINT *,'init mdw=', mdw
    177240
    178 !$OMP MASTER
    179 
    180 !    CALL gather(tr_seri, tr_seri_glo)
    181 !    IF (is_mpi_root) THEN
    182 !    IF (MAXVAL(tr_seri_glo).LT.1e-30) THEN
    183 
    184241    !--compute particle radius for a composition of 75% H2SO4 / 25% H2O at T=293K
    185242    DO bin_number=1, nbtr_bin
    … …  
    225282      &  lambda_int(NwvmaxSW+nwave_sw+nwave_lw+Nwv)
    226283
    227       OPEN (unit=11,file='h2so4_0.75_300.00_hummel_1988_p_q.dat')
    228 
    229284      IF (refr_ind_interpol) THEN
    230         DO nb=1,nb_lambda_h2so4
    231           READ(11,*) ref_ind(nb,:)
    232         ENDDO
     285
    233286        ilambda_max=ref_ind(1,2)/1.e6 !--in m
    234287        ilambda_min=ref_ind(nb_lambda_h2so4,2)/1.e6 !--in m
    … …  
    259312          ENDIF
    260313        ENDDO
    261       ELSE
    262         DO nb=1,nb_lambda_h2so4
    263           READ(11,*) wavenumber, ilambda, n_r_h2so4, n_i_h2so4
    264           ilambda=1000.*ilambda !wavelength in nm
    265           DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
    266             IF (ilambda/1.e9.GT.lambda_int(Nwv)) THEN !--- step function
    267               n_r(Nwv)=n_r_h2so4
    268               n_i(Nwv)=n_i_h2so4
     314
     315      ELSE  !-- no refr_ind_interpol, closest neighbour from upper wavelength
     316
     317        DO Nwv=1, NwvmaxSW+nwave_sw+nwave_lw+NwvmaxLW
     318          n_r(Nwv)=ref_ind(1,3)
     319          n_i(Nwv)=ref_ind(1,4)
     320          DO nb=2,nb_lambda_h2so4
     321            IF (ref_ind(nb,2)/1.e6.GT.lambda_int(Nwv)) THEN !--- step function
     322              n_r(Nwv)=ref_ind(nb,3)
     323              n_i(Nwv)=ref_ind(nb,4)
    269324            ENDIF 
    270325          ENDDO
    271326        ENDDO
    272327      ENDIF
    273       CLOSE(11)
    274328
    275329    !---Loop on wavelengths
    … …  
    400454
    401455      ENDDO  !--loop on wavelength
    402 
    403456
    404457    !---averaging over LMDZ wavebands
    … …  
    491544        piz_bin(nb,bin_number)=ss_w(nb)
    492545        cg_bin(nb,bin_number)=ss_g(nb)
    493 
    494546      ENDDO
    495547
    496548    ENDDO !loop over tracer bins
    497549
    498 !    ENDIF !mpi_root
    499 
    500 !$OMP END MASTER
    501   CALL bcast(alpha_bin)
    502   CALL bcast(piz_bin)
    503   CALL bcast(cg_bin)
    504 !$OMP BARRIER
    505 
    506 !    CALL scatter(alpha_bin, alpha_bin)
    507 !    CALL scatter(piz_bin, piz_bin)
    508 !    CALL scatter(cg_bin, cg_bin)
     550!!$OMP END MASTER
     551!  CALL bcast(alpha_bin)
     552!  CALL bcast(piz_bin)
     553!  CALL bcast(cg_bin)
     554!!$OMP BARRIER
    509555
    510556    !set to default values at first time step (tr_seri still zero)
    … …  
    515561    tau_strat_wave(:,:,:)=1.e-15
    516562
    517   ELSE
    518 
    519 !    CALL gather(tr_seri, tr_seri_glo)
     563  ELSE  !-- not debut
    520564
    521565  !--compute optical properties of actual size distribution (from tr_seri)

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

    r2700 r2704  
    7373    REAL,PARAMETER    :: xlat_sai=0.0           ! latitude of SAI in degree
    7474    REAL,PARAMETER    :: xlon_sai=120.35        ! longitude of SAI in degree
    75 
    76 !--be careful - this needs to be changed with resolution - here for 96x95 - seems to work for 48x36 as well
    77     REAL,PARAMETER    :: dlat=0.9474   ! d latitude in degree
    78     REAL,PARAMETER    :: dlon=1.875    ! d longitude in degree
    7975
    8076!--other local variables
    … …  
    9591    REAL                                   :: zdz                 ! thickness of atm. model layer in m
    9692    REAL,DIMENSION(klon,klev)              :: dens_aer            ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
     93    REAL                                   :: dlat, dlon          ! d latitude and d longitude of grid in degree
    9794
    9895    IF (is_mpi_root) THEN
    9996      PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
    10097    ENDIF
     98
     99    dlat=180./2./FLOAT(nbp_lat)   ! d latitude in degree
     100    dlon=360./2./FLOAT(nbp_lon)   ! d longitude in degree
    101101
    102102    DO it=1, nbtr_bin
    … …  
    168168        DO i=1,klon
    169169          !Pinatubo eruption at 15.14N, 120.35E
    170           IF  ( xlat(i).GT.xlat_vol-dlat .AND. xlat(i).LT.xlat_vol+dlat .AND. &
    171                 xlon(i).GT.xlon_vol-dlon .AND. xlon(i).LT.xlon_vol+dlon ) THEN
     170          IF  ( xlat(i).GE.xlat_vol-dlat .AND. xlat(i).LT.xlat_vol+dlat .AND. &
     171                xlon(i).GE.xlon_vol-dlon .AND. xlon(i).LT.xlon_vol+dlon ) THEN
    172172!         compute altLMDz
    173173            altLMDz(:)=0.0
    … …  
    209209!       IF  ((mth_cur==4 .AND. &
    210210!       SAI continuous emission o
    211         IF  ( xlat(i).GT.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. &
    212           &   xlon(i).GT.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
     211        IF  ( xlat(i).GE.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. &
     212          &   xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
    213213!         compute altLMDz
    214214          altLMDz(:)=0.0
Note: See TracChangeset for help on using the changeset viewer.