SUBROUTINE read_dust(debutphy, step, nbjour, dust_ec)
  USE dimphy
  USE lmdz_grid_phy
  USE lmdz_phys_para
  USE netcdf, ONLY: nf90_get_var, nf90_nowrite, nf90_open, nf90_inq_varid
  IMPLICIT NONE

  INCLUDE "dimensions.h"
  INCLUDE "paramet.h"

  INTEGER :: step, nbjour
  LOGICAL :: debutphy
  real :: dust_ec(klon)
  real :: dust_ec_glo(klon_glo)

  ! as      real dust_nc(iip1,jjp1)
  real :: dust_nc_glo(nbp_lon + 1, nbp_lat)
  integer :: rcode
  integer :: ncid1, varid1, ncid2, varid2

  save ncid1, varid1, ncid2, varid2
  !$OMP THREADPRIVATE(ncid1, varid1, ncid2, varid2)
  integer :: start(4), count(4), status
  integer :: i, j, ig

  !$OMP MASTER
  IF (is_mpi_root .AND. is_omp_root) THEN
    if (debutphy) then
      ncid1 = nf90_open('dust.nc', nf90_nowrite, rcode)
      varid1 = nf90_inq_varid(ncid1, 'EMISSION', rcode)
    endif

    start(1) = 1
    start(2) = 1
    start(4) = 0

    count(1) = nbp_lon + 1
    count(2) = nbp_lat
    count(3) = 1
    count(4) = 0

    start(3) = step

    status = nf90_get_var(ncid1, varid1, dust_nc_glo, start, count)

    ! CALL correctbid(iim,jjp1,dust_nc)
    CALL correctbid(nbp_lon, nbp_lat, dust_nc_glo)

    !--upside down + physical grid

    !--OB=change jjp1 to 1 here ;
    !----AS: OB's change is needed if lats N to S (90 to -90) in dust.nc
    !  dust_ec(1)=MAX(dust_nc(1,jjp1),0.0)
    dust_ec_glo(1) = MAX(dust_nc_glo(1, nbp_lat), 0.0)
    ig = 2
    ! DO j=2,jjm
    DO j = 2, nbp_lat - 1
      ! DO i = 1, iim
      DO i = 1, nbp_lon
        !--OB=change jjp1+1-j to j here
        ! dust_ec(ig)=MAX(dust_nc(i,jjp1+1-j),0.0)
        dust_ec_glo(ig) = MAX(dust_nc_glo(i, nbp_lat + 1 - j), 0.0)
        ig = ig + 1
      ENDDO
    ENDDO
    !--OB=change second 1 to jjp1 here
    dust_ec_glo(ig) = MAX(dust_nc_glo(1, 1), 0.0)
    ! end if master
  ENDIF
  !$OMP END MASTER
  !$OMP BARRIER
  CALL scatter(dust_ec_glo, dust_ec)


END SUBROUTINE read_dust