Ignore:
Timestamp:
Jul 24, 2024, 12:17:33 PM (4 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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/phylmd/readaerosol_mod.F90

    r5110 r5111  
    33MODULE readaerosol_mod
    44
    5   USE netcdf, ONLY: nf90_strerror,nf90_noerr,nf90_get_var,nf90_inq_varid,&
    6           nf90_inquire_dimension,nf90_inq_dimid,nf90_open,nf90_nowrite,nf90_close
    7 
    8   REAL, SAVE :: not_valid=-333.
    9  
     5  USE netcdf, ONLY: nf90_strerror, nf90_noerr, nf90_get_var, nf90_inq_varid, &
     6          nf90_inquire_dimension, nf90_inq_dimid, nf90_open, nf90_nowrite, nf90_close
     7  USE lmdz_abort_physic, ONLY: abort_physic
     8
     9  REAL, SAVE :: not_valid = -333.
     10
    1011  INTEGER, SAVE :: nbp_lon_src
    11 !$OMP THREADPRIVATE(nbp_lon_src) 
     12  !$OMP THREADPRIVATE(nbp_lon_src)
    1213  INTEGER, SAVE :: nbp_lat_src
    13 !$OMP THREADPRIVATE(nbp_lat_src) 
    14   REAL, ALLOCATABLE, SAVE    :: psurf_interp(:,:)
     14  !$OMP THREADPRIVATE(nbp_lat_src)
     15  REAL, ALLOCATABLE, SAVE :: psurf_interp(:, :)
    1516
    1617CONTAINS
    1718
    18 SUBROUTINE readaerosol(name_aero, type, filename, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
    19 
    20 !****************************************************************************************
    21 ! This routine will read the aersosol from file.
    22 
    23 ! Read a year data with get_aero_fromfile depending on aer_type :
    24 ! - actuel   : read year 1980
    25 ! - preind   : read natural data
    26 ! - scenario : read one or two years and do eventually linare time interpolation
    27 
    28 ! Return pointer, pt_out, to the year read or result from interpolation
    29 !****************************************************************************************
    30   USE dimphy
    31   USE print_control_mod, ONLY: lunout
    32 
    33   IMPLICIT NONE
    34 
    35   ! Input arguments
    36   CHARACTER(len=7), INTENT(IN) :: name_aero
    37   CHARACTER(len=*), INTENT(IN) :: type  ! actuel, annuel, scenario or preind
    38   CHARACTER(len=8), INTENT(IN) :: filename
    39   INTEGER, INTENT(IN)          :: iyr_in
    40 
    41   ! Output
    42   INTEGER, INTENT(OUT)            :: klev_src
    43   REAL, POINTER, DIMENSION(:)     :: pt_ap        ! Pointer for describing the vertical levels     
    44   REAL, POINTER, DIMENSION(:)     :: pt_b         ! Pointer for describing the vertical levels     
    45   REAL, POINTER, DIMENSION(:,:,:) :: pt_out       ! The massvar distributions, DIMENSION(klon, klev_src, 12)
    46   REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf  ! Surface pression for 12 months
    47   REAL, DIMENSION(klon,12), INTENT(OUT) :: load   ! Aerosol mass load in each column for 12 months
    48 
    49   ! Local variables
    50   CHARACTER(len=4)                :: cyear
    51   REAL, POINTER, DIMENSION(:,:,:) :: pt_2
    52   REAL, DIMENSION(klon,12)        :: psurf2, load2
    53   INTEGER                         :: iyr1, iyr2, klev_src2
    54   INTEGER                         :: it, k, i
    55   LOGICAL, PARAMETER              :: lonlyone=.FALSE.
    56 
    57 !****************************************************************************************
    58 ! Read data depending on aer_type
    59 
    60 !****************************************************************************************
    61 
    62   IF (type == 'actuel') THEN
    63 ! Read and return data for year 1980
    64 !****************************************************************************************
    65      cyear='1980'
    66      ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
    67      ! pt_out has dimensions (klon, klev_src, 12)
    68      CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
    69      
    70 
    71   ELSE IF (type == 'preind') THEN
    72 ! Read and return data from file with suffix .nat
    73 !****************************************************************************************     
    74      cyear='.nat'
    75      ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
    76      ! pt_out has dimensions (klon, klev_src, 12)
    77      CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
    78      
    79   ELSE IF (type == 'annuel') THEN
    80 ! Read and return data from scenario annual files
    81 !****************************************************************************************     
    82      WRITE(cyear,'(I4)') iyr_in
    83      WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,'   ',cyear
    84      ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month
    85      ! pt_out has dimensions (klon, klev_src, 12)
    86      CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
    87      
    88   ELSE IF (type == 'scenario') THEN
    89 ! Read data depending on actual year and interpolate if necessary
    90 !****************************************************************************************
    91      IF (iyr_in < 1850) THEN
    92         cyear='.nat'
    93         WRITE(lunout,*) 'get_aero 1 iyr_in=', iyr_in,'   ',cyear
     19  SUBROUTINE readaerosol(name_aero, type, filename, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
     20
     21    !****************************************************************************************
     22    ! This routine will read the aersosol from file.
     23
     24    ! Read a year data with get_aero_fromfile depending on aer_type :
     25    ! - actuel   : read year 1980
     26    ! - preind   : read natural data
     27    ! - scenario : read one or two years and do eventually linare time interpolation
     28
     29    ! Return pointer, pt_out, to the year read or result from interpolation
     30    !****************************************************************************************
     31    USE dimphy
     32    USE print_control_mod, ONLY: lunout
     33
     34    IMPLICIT NONE
     35
     36    ! Input arguments
     37    CHARACTER(len = 7), INTENT(IN) :: name_aero
     38    CHARACTER(len = *), INTENT(IN) :: type  ! actuel, annuel, scenario or preind
     39    CHARACTER(len = 8), INTENT(IN) :: filename
     40    INTEGER, INTENT(IN) :: iyr_in
     41
     42    ! Output
     43    INTEGER, INTENT(OUT) :: klev_src
     44    REAL, POINTER, DIMENSION(:) :: pt_ap        ! Pointer for describing the vertical levels
     45    REAL, POINTER, DIMENSION(:) :: pt_b         ! Pointer for describing the vertical levels
     46    REAL, POINTER, DIMENSION(:, :, :) :: pt_out       ! The massvar distributions, DIMENSION(klon, klev_src, 12)
     47    REAL, DIMENSION(klon, 12), INTENT(OUT) :: psurf  ! Surface pression for 12 months
     48    REAL, DIMENSION(klon, 12), INTENT(OUT) :: load   ! Aerosol mass load in each column for 12 months
     49
     50    ! Local variables
     51    CHARACTER(len = 4) :: cyear
     52    REAL, POINTER, DIMENSION(:, :, :) :: pt_2
     53    REAL, DIMENSION(klon, 12) :: psurf2, load2
     54    INTEGER :: iyr1, iyr2, klev_src2
     55    INTEGER :: it, k, i
     56    LOGICAL, PARAMETER :: lonlyone = .FALSE.
     57
     58    !****************************************************************************************
     59    ! Read data depending on aer_type
     60
     61    !****************************************************************************************
     62
     63    IF (type == 'actuel') THEN
     64      ! Read and return data for year 1980
     65      !****************************************************************************************
     66      cyear = '1980'
     67      ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
     68      ! pt_out has dimensions (klon, klev_src, 12)
     69      CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
     70
     71    ELSE IF (type == 'preind') THEN
     72      ! Read and return data from file with suffix .nat
     73      !****************************************************************************************
     74      cyear = '.nat'
     75      ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
     76      ! pt_out has dimensions (klon, klev_src, 12)
     77      CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
     78
     79    ELSE IF (type == 'annuel') THEN
     80      ! Read and return data from scenario annual files
     81      !****************************************************************************************
     82      WRITE(cyear, '(I4)') iyr_in
     83      WRITE(lunout, *) 'get_aero 3 iyr_in=', iyr_in, '   ', cyear
     84      ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month
     85      ! pt_out has dimensions (klon, klev_src, 12)
     86      CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
     87
     88    ELSE IF (type == 'scenario') THEN
     89      ! Read data depending on actual year and interpolate if necessary
     90      !****************************************************************************************
     91      IF (iyr_in < 1850) THEN
     92        cyear = '.nat'
     93        WRITE(lunout, *) 'get_aero 1 iyr_in=', iyr_in, '   ', cyear
    9494        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
    9595        ! pt_out has dimensions (klon, klev_src, 12)
    9696        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
    97        
    98      ELSE IF (iyr_in >= 2100) THEN
    99         cyear='2100'
    100         WRITE(lunout,*) 'get_aero 2 iyr_in=', iyr_in,'   ',cyear
     97
     98      ELSE IF (iyr_in >= 2100) THEN
     99        cyear = '2100'
     100        WRITE(lunout, *) 'get_aero 2 iyr_in=', iyr_in, '   ', cyear
    101101        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
    102102        ! pt_out has dimensions (klon, klev_src, 12)
    103103        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
    104        
    105      ELSE
     104
     105      ELSE
    106106        ! Read data from 2 decades and interpolate to actual year
    107107        ! a) from actual 10-yr-period
    108108        IF (iyr_in<1900) THEN
    109            iyr1 = 1850
    110            iyr2 = 1900
     109          iyr1 = 1850
     110          iyr2 = 1900
    111111        ELSE IF (iyr_in>=1900.AND.iyr_in<1920) THEN
    112            iyr1 = 1900
    113            iyr2 = 1920
    114         ELSE 
    115            iyr1 = INT(iyr_in/10)*10
    116            iyr2 = INT(1+iyr_in/10)*10
     112          iyr1 = 1900
     113          iyr2 = 1920
     114        ELSE
     115          iyr1 = INT(iyr_in / 10) * 10
     116          iyr2 = INT(1 + iyr_in / 10) * 10
    117117        ENDIF
    118        
    119         WRITE(cyear,'(I4)') iyr1
    120         WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,'   ',cyear
     118
     119        WRITE(cyear, '(I4)') iyr1
     120        WRITE(lunout, *) 'get_aero 3 iyr_in=', iyr_in, '   ', cyear
    121121        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
    122122        ! pt_out has dimensions (klon, klev_src, 12)
    123123        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
    124        
     124
    125125        ! If to read two decades:
    126         IF (.NOT.lonlyone) THEN
    127            
    128            ! b) from the next following one
    129            WRITE(cyear,'(I4)') iyr2
    130            WRITE(lunout,*) 'get_aero 4 iyr_in=', iyr_in,'   ',cyear
    131            
    132            NULLIFY(pt_2)
    133            ! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month
    134            ! pt_2 has dimensions (klon, klev_src, 12)
    135            CALL get_aero_fromfile(name_aero, cyear, filename, klev_src2, pt_ap, pt_b, pt_2, psurf2, load2)
    136            ! Test for same number of vertical levels
    137            IF (klev_src /= klev_src2) THEN
    138               WRITE(lunout,*) 'Two aerosols files with different number of vertical levels is not allowded'
    139               CALL abort_physic('readaersosol','Error in number of vertical levels',1)
    140            END IF
    141            
    142            ! Linare interpolate to the actual year:
    143            DO it=1,12
    144               DO k=1,klev_src
    145                  DO i = 1, klon
    146                     pt_out(i,k,it) = &
    147                          pt_out(i,k,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * &
    148                          (pt_out(i,k,it) - pt_2(i,k,it))
    149                  END DO
     126        IF (.NOT.lonlyone) THEN
     127
     128          ! b) from the next following one
     129          WRITE(cyear, '(I4)') iyr2
     130          WRITE(lunout, *) 'get_aero 4 iyr_in=', iyr_in, '   ', cyear
     131
     132          NULLIFY(pt_2)
     133          ! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month
     134          ! pt_2 has dimensions (klon, klev_src, 12)
     135          CALL get_aero_fromfile(name_aero, cyear, filename, klev_src2, pt_ap, pt_b, pt_2, psurf2, load2)
     136          ! Test for same number of vertical levels
     137          IF (klev_src /= klev_src2) THEN
     138            WRITE(lunout, *) 'Two aerosols files with different number of vertical levels is not allowded'
     139            CALL abort_physic('readaersosol', 'Error in number of vertical levels', 1)
     140          END IF
     141
     142          ! Linare interpolate to the actual year:
     143          DO it = 1, 12
     144            DO k = 1, klev_src
     145              DO i = 1, klon
     146                pt_out(i, k, it) = &
     147                        pt_out(i, k, it) - REAL(iyr_in - iyr1) / REAL(iyr2 - iyr1) * &
     148                                (pt_out(i, k, it) - pt_2(i, k, it))
    150149              END DO
    151 
    152               DO i = 1, klon
    153                  psurf(i,it) = &
    154                       psurf(i,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * &
    155                       (psurf(i,it) - psurf2(i,it))
    156 
    157                  load(i,it) = &
    158                       load(i,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * &
    159                       (load(i,it) - load2(i,it))
    160               END DO
    161            END DO
    162 
    163            ! Deallocate pt_2 no more needed
    164            DEALLOCATE(pt_2)
    165            
     150            END DO
     151
     152            DO i = 1, klon
     153              psurf(i, it) = &
     154                      psurf(i, it) - REAL(iyr_in - iyr1) / REAL(iyr2 - iyr1) * &
     155                              (psurf(i, it) - psurf2(i, it))
     156
     157              load(i, it) = &
     158                      load(i, it) - REAL(iyr_in - iyr1) / REAL(iyr2 - iyr1) * &
     159                              (load(i, it) - load2(i, it))
     160            END DO
     161          END DO
     162
     163          ! Deallocate pt_2 no more needed
     164          DEALLOCATE(pt_2)
     165
    166166        END IF ! lonlyone
    167      END IF ! iyr_in .LT. 1850
    168 
    169   ELSE
    170      WRITE(lunout,*)'This option is not implemented : aer_type = ', type,' name_aero=',name_aero
    171      CALL abort_physic('readaerosol','Error : aer_type parameter not accepted',1)
    172   END IF ! type
    173 
    174 
    175 END SUBROUTINE readaerosol
    176 
    177 
    178 SUBROUTINE init_aero_fromfile(flag_aerosol, aerosol_couple)
    179   USE lmdz_phys_para
    180   USE lmdz_grid_phy, ONLY: grid_type, unstructured
    181   USE lmdz_xios
    182   IMPLICIT NONE
    183 
    184   INTEGER, INTENT(IN) :: flag_aerosol
    185   LOGICAL, INTENT(IN) :: aerosol_couple
    186  
    187   REAL,ALLOCATABLE :: lat_src(:)
    188   REAL,ALLOCATABLE :: lon_src(:)
    189   CHARACTER(LEN=*),PARAMETER :: file_aerosol='aerosols.nat.nc'
    190   CHARACTER(LEN=*),PARAMETER :: file_so4='so4.nat.nc'
    191   INTEGER :: klev_src
    192   INTEGER :: ierr ,ncid, dimID, varid
    193   REAL :: null_array(0)
    194 
    195   IF (using_xios) THEN
    196     IF (flag_aerosol>0 .AND. grid_type==unstructured .AND. (.NOT. aerosol_couple) ) THEN
    197  
    198       IF (is_omp_root) THEN
    199  
    200         IF (is_mpi_root) THEN
    201    
    202           IF (nf90_open(TRIM(file_aerosol), nf90_nowrite, ncid) /= nf90_noerr) THEN
    203             CALL check_err( nf90_open(TRIM(file_so4), nf90_nowrite, ncid), "pb open "//trim(file_so4) )
     167      END IF ! iyr_in .LT. 1850
     168
     169    ELSE
     170      WRITE(lunout, *)'This option is not implemented : aer_type = ', type, ' name_aero=', name_aero
     171      CALL abort_physic('readaerosol', 'Error : aer_type parameter not accepted', 1)
     172    END IF ! type
     173
     174  END SUBROUTINE readaerosol
     175
     176
     177  SUBROUTINE init_aero_fromfile(flag_aerosol, aerosol_couple)
     178    USE lmdz_phys_para
     179    USE lmdz_grid_phy, ONLY: grid_type, unstructured
     180    USE lmdz_xios
     181    IMPLICIT NONE
     182
     183    INTEGER, INTENT(IN) :: flag_aerosol
     184    LOGICAL, INTENT(IN) :: aerosol_couple
     185
     186    REAL, ALLOCATABLE :: lat_src(:)
     187    REAL, ALLOCATABLE :: lon_src(:)
     188    CHARACTER(LEN = *), PARAMETER :: file_aerosol = 'aerosols.nat.nc'
     189    CHARACTER(LEN = *), PARAMETER :: file_so4 = 'so4.nat.nc'
     190    INTEGER :: klev_src
     191    INTEGER :: ierr, ncid, dimID, varid
     192    REAL :: null_array(0)
     193
     194    IF (using_xios) THEN
     195      IF (flag_aerosol>0 .AND. grid_type==unstructured .AND. (.NOT. aerosol_couple)) THEN
     196
     197        IF (is_omp_root) THEN
     198
     199          IF (is_mpi_root) THEN
     200
     201            IF (nf90_open(TRIM(file_aerosol), nf90_nowrite, ncid) /= nf90_noerr) THEN
     202              CALL check_err(nf90_open(TRIM(file_so4), nf90_nowrite, ncid), "pb open " // trim(file_so4))
     203            ENDIF
     204
     205            ! Read and test longitudes
     206            CALL check_err(nf90_inq_dimid(ncid, "lon", dimID), "pb inq dim lon")
     207            CALL check_err(nf90_inquire_dimension(ncid, dimID, len = nbp_lon_src), "pb inq dim lon")
     208            CALL check_err(nf90_inq_varid(ncid, 'lon', varid), "pb inq lon")
     209            ALLOCATE(lon_src(nbp_lon_src))
     210            CALL check_err(nf90_get_var(ncid, varid, lon_src(:)), "pb get lon")
     211
     212            ! Read and test latitudes
     213            CALL check_err(nf90_inq_dimid(ncid, "lat", dimID), "pb inq dim lat")
     214            CALL check_err(nf90_inquire_dimension(ncid, dimID, len = nbp_lat_src), "pb inq dim lat")
     215            CALL check_err(nf90_inq_varid(ncid, 'lat', varid), "pb inq lat")
     216            ALLOCATE(lat_src(nbp_lat_src))
     217            CALL check_err(nf90_get_var(ncid, varid, lat_src(:)), "pb get lat")
     218            IF (nf90_inq_dimid(ncid, 'lev', dimid) /= nf90_noerr) THEN
     219              IF (nf90_inq_dimid(ncid, 'presnivs', dimid)/= nf90_noerr) THEN
     220                CALL check_err(nf90_inq_dimid(ncid, 'PRESNIVS', dimid), 'dimension lev,PRESNIVS or presnivs not in file')
     221              ENDIF
     222            ENDIF
     223            CALL check_err(nf90_inquire_dimension(ncid, dimid, len = klev_src), "pb inq dim for PRESNIVS or lev")
     224            CALL check_err(nf90_close(ncid), "pb in close")
    204225          ENDIF
    205226
    206           ! Read and test longitudes
    207           CALL check_err( nf90_inq_dimid(ncid, "lon", dimID),"pb inq dim lon")
    208           CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lon_src),"pb inq dim lon")
    209           CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
    210           ALLOCATE(lon_src(nbp_lon_src))
    211           CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
    212 
    213           ! Read and test latitudes
    214           CALL check_err( nf90_inq_dimid(ncid, "lat", dimID),"pb inq dim lat")
    215           CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lat_src),"pb inq dim lat")
    216           CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
    217           ALLOCATE(lat_src(nbp_lat_src))
    218           CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
    219           IF (nf90_inq_dimid(ncid, 'lev', dimid) /= nf90_noerr) THEN
    220             IF (nf90_inq_dimid(ncid, 'presnivs', dimid)/= nf90_noerr) THEN
    221                CALL check_err(nf90_inq_dimid(ncid, 'PRESNIVS', dimid),'dimension lev,PRESNIVS or presnivs not in file')
    222             ENDIF
     227          CALL bcast_mpi(nbp_lat_src)
     228          CALL bcast_mpi(nbp_lon_src)
     229          CALL bcast_mpi(klev_src)
     230
     231          IF (is_mpi_root) THEN
     232            CALL xios_set_domain_attr("domain_aerosol", nj_glo = nbp_lat_src, nj = nbp_lat_src, jbegin = 0, latvalue_1d = lat_src)
     233            CALL xios_set_domain_attr("domain_aerosol", ni_glo = nbp_lon_src, ni = nbp_lon_src, ibegin = 0, lonvalue_1d = lon_src)
     234          ELSE
     235            CALL xios_set_domain_attr("domain_aerosol", nj_glo = nbp_lat_src, nj = 0, jbegin = 0, latvalue_1d = null_array)
     236            CALL xios_set_domain_attr("domain_aerosol", ni_glo = nbp_lon_src, ni = 0, ibegin = 0, lonvalue_1d = null_array)
    223237          ENDIF
    224           CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
    225           CALL check_err( nf90_close(ncid),"pb in close" )   
     238          CALL xios_set_axis_attr("axis_aerosol", n_glo = klev_src)
     239          CALL xios_set_fieldgroup_attr("aerosols", enabled = .TRUE.)
     240
    226241        ENDIF
    227242
    228         CALL bcast_mpi(nbp_lat_src)
    229         CALL bcast_mpi(nbp_lon_src)
    230         CALL bcast_mpi(klev_src)
    231 
    232         IF (is_mpi_root ) THEN
    233           CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=nbp_lat_src, jbegin=0, latvalue_1d=lat_src)
    234           CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=nbp_lon_src, ibegin=0, lonvalue_1d=lon_src)
    235         ELSE
    236           CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=0, jbegin=0, latvalue_1d=null_array )
    237           CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=0, ibegin=0, lonvalue_1d=null_array)
    238         ENDIF
    239         CALL xios_set_axis_attr("axis_aerosol",n_glo=klev_src)
    240         CALL xios_set_fieldgroup_attr("aerosols", enabled=.TRUE.)
    241  
    242243      ENDIF
    243    
    244     ENDIF   
    245   ENDIF !using_xios
    246 END SUBROUTINE init_aero_fromfile
    247 
    248 
    249    
     244    ENDIF !using_xios
     245  END SUBROUTINE init_aero_fromfile
     246
    250247
    251248  SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, pt_year, psurf_out, load_out)
    252 !****************************************************************************************
    253 ! Read 12 month aerosol from file and distribute to local process on physical grid.
    254 ! Vertical levels, klev_src, may differ from model levels if new file format.
    255 
    256 ! For mpi_root and master thread :
    257 ! 1) Open file
    258 ! 2) Find vertical dimension klev_src
    259 ! 3) Read field month by month
    260 ! 4) Close file 
    261 ! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
    262 !     - Also the levels and the latitudes have to be inversed
    263 
    264 ! For all processes and threads :
    265 ! 6) Scatter global field(klon_glo) to local process domain(klon)
    266 ! 7) Test for negative values
    267 !****************************************************************************************
     249    !****************************************************************************************
     250    ! Read 12 month aerosol from file and distribute to local process on physical grid.
     251    ! Vertical levels, klev_src, may differ from model levels if new file format.
     252
     253    ! For mpi_root and master thread :
     254    ! 1) Open file
     255    ! 2) Find vertical dimension klev_src
     256    ! 3) Read field month by month
     257    ! 4) Close file
     258    ! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
     259    !     - Also the levels and the latitudes have to be inversed
     260
     261    ! For all processes and threads :
     262    ! 6) Scatter global field(klon_glo) to local process domain(klon)
     263    ! 7) Test for negative values
     264    !****************************************************************************************
    268265
    269266    USE dimphy
    270     USE lmdz_grid_phy, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, klon_glo, &
    271                                  grid2Dto1D_glo, grid_type, unstructured
     267    USE lmdz_grid_phy, ONLY: nbp_lon_ => nbp_lon, nbp_lat_ => nbp_lat, klon_glo, &
     268            grid2Dto1D_glo, grid_type, unstructured
    272269    USE lmdz_phys_para
    273270    USE iophy, ONLY: io_lon, io_lat
    … …  
    275272    USE lmdz_xios
    276273    IMPLICIT NONE
    277      
    278 ! Input argumets
    279     CHARACTER(len=7), INTENT(IN)          :: varname
    280     CHARACTER(len=4), INTENT(IN)          :: cyr
    281     CHARACTER(len=8), INTENT(IN)          :: filename
    282 
    283 ! Output arguments
    284     INTEGER, INTENT(OUT)                  :: klev_src     ! Number of vertical levels in file
    285     REAL, POINTER, DIMENSION(:)           :: pt_ap        ! Pointer for describing the vertical levels     
    286     REAL, POINTER, DIMENSION(:)           :: pt_b         ! Pointer for describing the vertical levels     
    287     REAL, POINTER, DIMENSION(:,:,:)      :: pt_year      ! Pointer-variabale from file, 12 month, grid : klon,klev_src
    288     REAL, POINTER, DIMENSION(:,:,:)      :: pt_year_mpi  ! Pointer-variabale from file, 12 month, grid : klon,klev_src
    289     REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out    ! Surface pression for 12 months
    290     REAL, DIMENSION(klon_mpi,12)          :: psurf_out_mpi    ! Surface pression for 12 months
    291     REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out     ! Aerosol mass load in each column
    292     REAL, DIMENSION(klon_mpi,12)          :: load_out_mpi     ! Aerosol mass load in each column
    293     INTEGER                               :: nbr_tsteps   ! number of month in file read
    294 
    295 ! Local variables
    296     CHARACTER(len=30)    :: fname
    297     CHARACTER(len=30)    :: cvar
    298     INTEGER               :: ncid, dimid, varid
    299     INTEGER               :: imth, i, j, k, ierr
    300     REAL                  :: npole, spole
    301     REAL, ALLOCATABLE, DIMENSION(:,:,:)  :: varmth
    302     REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear       ! Global variable read from file, 12 month
    303     REAL, ALLOCATABLE, DIMENSION(:,:,:)  :: varyear_glo1D !(klon_glo, klev_src, 12)
    304     REAL, ALLOCATABLE, DIMENSION(:)       :: varktmp
    305 
    306     REAL,  ALLOCATABLE                    :: psurf_glo2D(:,:,:)   ! Surface pression for 12 months on dynamics global grid
    307     REAL, DIMENSION(klon_glo,12)          :: psurf_glo1D   ! -"- on physical global grid
    308     REAL,  ALLOCATABLE                    :: load_glo2D(:,:,:)    ! Load for 12 months on dynamics global grid
    309     REAL, DIMENSION(klon_glo,12)          :: load_glo1D    ! -"- on physical global grid
    310     REAL, ALLOCATABLE, DIMENSION(:,:)    :: vartmp
    311     REAL, ALLOCATABLE,DIMENSION(:)        :: lon_src              ! longitudes in file
    312     REAL, ALLOCATABLE,DIMENSION(:)        :: lat_src, lat_src_inv ! latitudes in file
    313     LOGICAL                               :: new_file             ! true if new file format detected
    314     LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
    315     INTEGER                               :: nbp_lon, nbp_lat
    316     LOGICAL,SAVE                          :: first=.TRUE.
    317 !$OMP THREADPRIVATE(first)
    318    
     274
     275    ! Input argumets
     276    CHARACTER(len = 7), INTENT(IN) :: varname
     277    CHARACTER(len = 4), INTENT(IN) :: cyr
     278    CHARACTER(len = 8), INTENT(IN) :: filename
     279
     280    ! Output arguments
     281    INTEGER, INTENT(OUT) :: klev_src     ! Number of vertical levels in file
     282    REAL, POINTER, DIMENSION(:) :: pt_ap        ! Pointer for describing the vertical levels
     283    REAL, POINTER, DIMENSION(:) :: pt_b         ! Pointer for describing the vertical levels
     284    REAL, POINTER, DIMENSION(:, :, :) :: pt_year      ! Pointer-variabale from file, 12 month, grid : klon,klev_src
     285    REAL, POINTER, DIMENSION(:, :, :) :: pt_year_mpi  ! Pointer-variabale from file, 12 month, grid : klon,klev_src
     286    REAL, DIMENSION(klon, 12), INTENT(OUT) :: psurf_out    ! Surface pression for 12 months
     287    REAL, DIMENSION(klon_mpi, 12) :: psurf_out_mpi    ! Surface pression for 12 months
     288    REAL, DIMENSION(klon, 12), INTENT(OUT) :: load_out     ! Aerosol mass load in each column
     289    REAL, DIMENSION(klon_mpi, 12) :: load_out_mpi     ! Aerosol mass load in each column
     290    INTEGER :: nbr_tsteps   ! number of month in file read
     291
     292    ! Local variables
     293    CHARACTER(len = 30) :: fname
     294    CHARACTER(len = 30) :: cvar
     295    INTEGER :: ncid, dimid, varid
     296    INTEGER :: imth, i, j, k, ierr
     297    REAL :: npole, spole
     298    REAL, ALLOCATABLE, DIMENSION(:, :, :) :: varmth
     299    REAL, ALLOCATABLE, DIMENSION(:, :, :, :) :: varyear       ! Global variable read from file, 12 month
     300    REAL, ALLOCATABLE, DIMENSION(:, :, :) :: varyear_glo1D !(klon_glo, klev_src, 12)
     301    REAL, ALLOCATABLE, DIMENSION(:) :: varktmp
     302
     303    REAL, ALLOCATABLE :: psurf_glo2D(:, :, :)   ! Surface pression for 12 months on dynamics global grid
     304    REAL, DIMENSION(klon_glo, 12) :: psurf_glo1D   ! -"- on physical global grid
     305    REAL, ALLOCATABLE :: load_glo2D(:, :, :)    ! Load for 12 months on dynamics global grid
     306    REAL, DIMENSION(klon_glo, 12) :: load_glo1D    ! -"- on physical global grid
     307    REAL, ALLOCATABLE, DIMENSION(:, :) :: vartmp
     308    REAL, ALLOCATABLE, DIMENSION(:) :: lon_src              ! longitudes in file
     309    REAL, ALLOCATABLE, DIMENSION(:) :: lat_src, lat_src_inv ! latitudes in file
     310    LOGICAL :: new_file             ! true if new file format detected
     311    LOGICAL :: invert_lat           ! true if the field has to be inverted for latitudes
     312    INTEGER :: nbp_lon, nbp_lat
     313    LOGICAL, SAVE :: first = .TRUE.
     314    !$OMP THREADPRIVATE(first)
     315
    319316    IF (grid_type==unstructured) THEN
    320       nbp_lon=nbp_lon_src
    321       nbp_lat=nbp_lat_src
     317      nbp_lon = nbp_lon_src
     318      nbp_lat = nbp_lat_src
    322319    ELSE
    323       nbp_lon=nbp_lon_
    324       nbp_lat=nbp_lat_
     320      nbp_lon = nbp_lon_
     321      nbp_lat = nbp_lat_
    325322    ENDIF
    326    
     323
    327324    IF (is_mpi_root) THEN
    328    
    329       ALLOCATE(psurf_glo2D(nbp_lon,nbp_lat,12))
    330       ALLOCATE(load_glo2D(nbp_lon,nbp_lat,12))
    331       ALLOCATE(vartmp(nbp_lon,nbp_lat))
     325
     326      ALLOCATE(psurf_glo2D(nbp_lon, nbp_lat, 12))
     327      ALLOCATE(load_glo2D(nbp_lon, nbp_lat, 12))
     328      ALLOCATE(vartmp(nbp_lon, nbp_lat))
    332329      ALLOCATE(lon_src(nbp_lon))
    333330      ALLOCATE(lat_src(nbp_lat))
    334331      ALLOCATE(lat_src_inv(nbp_lat))
    335332    ELSE
    336       ALLOCATE(varyear(0,0,0,0))
    337       ALLOCATE(psurf_glo2D(0,0,0))
    338       ALLOCATE(load_glo2D(0,0,0))
     333      ALLOCATE(varyear(0, 0, 0, 0))
     334      ALLOCATE(psurf_glo2D(0, 0, 0))
     335      ALLOCATE(load_glo2D(0, 0, 0))
    339336    ENDIF
    340            
     337
    341338    ! Deallocate pointers
    342339    IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
    … …  
    345342    IF (is_mpi_root .AND. is_omp_root) THEN
    346343
    347 ! 1) Open file
    348 !****************************************************************************************
    349 ! Add suffix to filename
    350        fname = trim(filename)//cyr//'.nc'
    351  
    352        WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname)
    353        CALL check_err( nf90_open(TRIM(fname), nf90_nowrite, ncid), "pb open "//trim(fname) )
    354 
    355 
    356        IF (grid_type/=unstructured) THEN
    357 
    358 ! Test for equal longitudes and latitudes in file and model
    359 !****************************************************************************************
    360        ! Read and test longitudes
    361          CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
    362          CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
    363        
    364          IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
    365             WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
    366             WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
    367             WRITE(lunout,*) 'longitudes in model :', io_lon
    368          
    369             CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
    370          END IF
    371 
    372          ! Read and test latitudes
    373          CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
    374          CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
    375 
    376          ! Invert source latitudes
    377          DO j = 1, nbp_lat
    378             lat_src_inv(j) = lat_src(nbp_lat +1 -j)
    379          END DO
    380 
    381          IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
    382             ! Latitudes are the same
    383             invert_lat=.FALSE.
    384          ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
    385             ! Inverted source latitudes correspond to model latitudes
    386             WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
    387             invert_lat=.TRUE.
    388          ELSE
    389             WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
    390             WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
    391             WRITE(lunout,*) 'latitudes in model :', io_lat
    392             CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
    393          END IF
    394        ENDIF
    395 
    396 ! 2) Check if old or new file is avalabale.
    397 !    New type of file should contain the dimension 'lev'
    398 !    Old type of file should contain the dimension 'PRESNIVS'
    399 !****************************************************************************************
    400        ierr = nf90_inq_dimid(ncid, 'lev', dimid)
    401        IF (ierr /= nf90_noerr) THEN
    402           ! Coordinate axe lev not found. Check for presnivs.
    403           ierr = nf90_inq_dimid(ncid, 'presnivs', dimid)
     344      ! 1) Open file
     345      !****************************************************************************************
     346      ! Add suffix to filename
     347      fname = trim(filename) // cyr // '.nc'
     348
     349      WRITE(lunout, *) 'reading variable ', TRIM(varname), ' in file ', TRIM(fname)
     350      CALL check_err(nf90_open(TRIM(fname), nf90_nowrite, ncid), "pb open " // trim(fname))
     351
     352      IF (grid_type/=unstructured) THEN
     353
     354        ! Test for equal longitudes and latitudes in file and model
     355        !****************************************************************************************
     356        ! Read and test longitudes
     357        CALL check_err(nf90_inq_varid(ncid, 'lon', varid), "pb inq lon")
     358        CALL check_err(nf90_get_var(ncid, varid, lon_src(:)), "pb get lon")
     359
     360        IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
     361          WRITE(lunout, *) 'Problem in longitudes read from file : ', TRIM(fname)
     362          WRITE(lunout, *) 'longitudes in file ', TRIM(fname), ' : ', lon_src
     363          WRITE(lunout, *) 'longitudes in model :', io_lon
     364
     365          CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model', 1)
     366        END IF
     367
     368        ! Read and test latitudes
     369        CALL check_err(nf90_inq_varid(ncid, 'lat', varid), "pb inq lat")
     370        CALL check_err(nf90_get_var(ncid, varid, lat_src(:)), "pb get lat")
     371
     372        ! Invert source latitudes
     373        DO j = 1, nbp_lat
     374          lat_src_inv(j) = lat_src(nbp_lat + 1 - j)
     375        END DO
     376
     377        IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
     378          ! Latitudes are the same
     379          invert_lat = .FALSE.
     380        ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
     381          ! Inverted source latitudes correspond to model latitudes
     382          WRITE(lunout, *) 'latitudes will be inverted for file : ', TRIM(fname)
     383          invert_lat = .TRUE.
     384        ELSE
     385          WRITE(lunout, *) 'Problem in latitudes read from file : ', TRIM(fname)
     386          WRITE(lunout, *) 'latitudes in file ', TRIM(fname), ' : ', lat_src
     387          WRITE(lunout, *) 'latitudes in model :', io_lat
     388          CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model', 1)
     389        END IF
     390      ENDIF
     391
     392      ! 2) Check if old or new file is avalabale.
     393      !    New type of file should contain the dimension 'lev'
     394      !    Old type of file should contain the dimension 'PRESNIVS'
     395      !****************************************************************************************
     396      ierr = nf90_inq_dimid(ncid, 'lev', dimid)
     397      IF (ierr /= nf90_noerr) THEN
     398        ! Coordinate axe lev not found. Check for presnivs.
     399        ierr = nf90_inq_dimid(ncid, 'presnivs', dimid)
     400        IF (ierr /= nf90_noerr) THEN
     401          ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid)
    404402          IF (ierr /= nf90_noerr) THEN
    405              ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid)
    406              IF (ierr /= nf90_noerr) THEN
    407                 ! Dimension PRESNIVS not found either
    408                 CALL abort_physic('get_aero_fromfile', 'dimension lev,PRESNIVS or presnivs not in file',1)
    409              ELSE
    410                 ! Old file found
    411                 new_file=.FALSE.
    412                 WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done'
    413              END IF
     403            ! Dimension PRESNIVS not found either
     404            CALL abort_physic('get_aero_fromfile', 'dimension lev,PRESNIVS or presnivs not in file', 1)
    414405          ELSE
    415              ! New file found
    416              new_file=.TRUE.
    417              WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
    418           ENDIF
    419        ELSE
     406            ! Old file found
     407            new_file = .FALSE.
     408            WRITE(lunout, *) 'Vertical interpolation for ', TRIM(varname), ' will not be done'
     409          END IF
     410        ELSE
    420411          ! New file found
    421           new_file=.TRUE.
    422           WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
    423        END IF
    424        
    425 ! 2) Find vertical dimension klev_src
    426 !****************************************************************************************
    427        CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
    428        
    429      ! Allocate variables depending on the number of vertical levels
    430        ALLOCATE(varmth(nbp_lon,nbp_lat, klev_src), varyear(nbp_lon,nbp_lat, klev_src, 12), stat=ierr)
    431        IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 1',1)
    432 
    433        ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr)
    434        IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 2',1)
    435 
    436 ! 3) Read all variables from file
    437 !    There is 2 options for the file structure :
    438 !    new_file=TRUE  : read varyear, ps, pt_ap and pt_b
    439 !    new_file=FALSE : read varyear month by month
    440 !****************************************************************************************
    441 
    442        IF (new_file) THEN
    443 ! ++) Check number of month in file opened
    444 !**************************************************************************************************
    445        ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
    446        if (ierr /= nf90_noerr) THEN
     412          new_file = .TRUE.
     413          WRITE(lunout, *) 'Vertical interpolation for ', TRIM(varname), ' will be done'
     414        ENDIF
     415      ELSE
     416        ! New file found
     417        new_file = .TRUE.
     418        WRITE(lunout, *) 'Vertical interpolation for ', TRIM(varname), ' will be done'
     419      END IF
     420
     421      ! 2) Find vertical dimension klev_src
     422      !****************************************************************************************
     423      CALL check_err(nf90_inquire_dimension(ncid, dimid, len = klev_src), "pb inq dim for PRESNIVS or lev")
     424
     425      ! Allocate variables depending on the number of vertical levels
     426      ALLOCATE(varmth(nbp_lon, nbp_lat, klev_src), varyear(nbp_lon, nbp_lat, klev_src, 12), stat = ierr)
     427      IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 1', 1)
     428
     429      ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat = ierr)
     430      IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 2', 1)
     431
     432      ! 3) Read all variables from file
     433      !    There is 2 options for the file structure :
     434      !    new_file=TRUE  : read varyear, ps, pt_ap and pt_b
     435      !    new_file=FALSE : read varyear month by month
     436      !****************************************************************************************
     437
     438      IF (new_file) THEN
     439        ! ++) Check number of month in file opened
     440        !**************************************************************************************************
     441        ierr = nf90_inq_dimid(ncid, 'TIME', dimid)
     442        if (ierr /= nf90_noerr) THEN
    447443          ierr = nf90_inq_dimid(ncid, 'time_counter', dimid)
    448        ENDIF
    449        CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps),"pb inq dim TIME or time_counter" )
    450 !       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
    451        IF (nbr_tsteps /= 12 ) THEN
    452     CALL abort_physic('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)' &
    453                      ,1)
    454        ENDIF
    455 
    456 ! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
    457 !****************************************************************************************
     444        ENDIF
     445        CALL check_err(nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps), "pb inq dim TIME or time_counter")
     446        !       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
     447        IF (nbr_tsteps /= 12) THEN
     448          CALL abort_physic('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)' &
     449                  , 1)
     450        ENDIF
     451
     452        ! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
     453        !****************************************************************************************
     454        ! Get variable id
     455        !CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) )
     456        print *, 'readaerosol ', TRIM(varname)
     457        IF (nf90_inq_varid(ncid, TRIM(varname), varid) /= nf90_noerr) THEN
     458          ! Variable is not there
     459          WRITE(lunout, *) 'Attention ' // TRIM(varname) // ' is not in aerosol input file'
     460          varyear(:, :, :, :) = 0.0
     461        ELSE
     462          ! Get the variable
     463          CALL check_err(nf90_get_var(ncid, varid, varyear(:, :, :, :)), "pb get var " // TRIM(varname))
     464        ENDIF
     465
     466        ! ++) Read surface pression, 12 month in one variable
     467        !****************************************************************************************
     468        ! Get variable id
     469        CALL check_err(nf90_inq_varid(ncid, "ps", varid), "pb inq var ps")
     470        ! Get the variable
     471        CALL check_err(nf90_get_var(ncid, varid, psurf_glo2D), "pb get var ps")
     472
     473        ! ++) Read mass load, 12 month in one variable
     474        !****************************************************************************************
     475        ! Get variable id
     476        !CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname))
     477        IF (nf90_inq_varid(ncid, "load_" // TRIM(varname), varid) /= nf90_noerr) THEN
     478          WRITE(lunout, *) 'Attention load_' // TRIM(varname) // ' is not in aerosol input file'
     479          load_glo2D(:, :, :) = 0.0
     480        ELSE
     481          ! Get the variable
     482          CALL check_err(nf90_get_var(ncid, varid, load_glo2D), "pb get var load_" // TRIM(varname))
     483        ENDIF
     484
     485        ! ++) Read ap
     486        !****************************************************************************************
     487        ! Get variable id
     488        CALL check_err(nf90_inq_varid(ncid, "ap", varid), "pb inq var ap")
     489        ! Get the variable
     490        CALL check_err(nf90_get_var(ncid, varid, pt_ap), "pb get var ap")
     491
     492        ! ++) Read b
     493        !****************************************************************************************
     494        ! Get variable id
     495        CALL check_err(nf90_inq_varid(ncid, "b", varid), "pb inq var b")
     496        ! Get the variable
     497        CALL check_err(nf90_get_var(ncid, varid, pt_b), "pb get var b")
     498
     499      ELSE  ! old file
     500
     501        ! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
     502        !****************************************************************************************
     503        DO imth = 1, 12
     504          IF (imth==1) THEN
     505            cvar = TRIM(varname) // 'JAN'
     506          ELSE IF (imth==2) THEN
     507            cvar = TRIM(varname) // 'FEB'
     508          ELSE IF (imth==3) THEN
     509            cvar = TRIM(varname) // 'MAR'
     510          ELSE IF (imth==4) THEN
     511            cvar = TRIM(varname) // 'APR'
     512          ELSE IF (imth==5) THEN
     513            cvar = TRIM(varname) // 'MAY'
     514          ELSE IF (imth==6) THEN
     515            cvar = TRIM(varname) // 'JUN'
     516          ELSE IF (imth==7) THEN
     517            cvar = TRIM(varname) // 'JUL'
     518          ELSE IF (imth==8) THEN
     519            cvar = TRIM(varname) // 'AUG'
     520          ELSE IF (imth==9) THEN
     521            cvar = TRIM(varname) // 'SEP'
     522          ELSE IF (imth==10) THEN
     523            cvar = TRIM(varname) // 'OCT'
     524          ELSE IF (imth==11) THEN
     525            cvar = TRIM(varname) // 'NOV'
     526          ELSE IF (imth==12) THEN
     527            cvar = TRIM(varname) // 'DEC'
     528          END IF
     529
    458530          ! Get variable id
    459           !CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) )
    460           print *,'readaerosol ', TRIM(varname)
    461           IF ( nf90_inq_varid(ncid, TRIM(varname), varid) /= nf90_noerr ) THEN
    462             ! Variable is not there
    463             WRITE(lunout,*) 'Attention '//TRIM(varname)//' is not in aerosol input file'
    464             varyear(:,:,:,:)=0.0
    465           ELSE
    466             ! Get the variable
    467             CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)),"pb get var "//TRIM(varname) )
    468           ENDIF
    469          
    470 ! ++) Read surface pression, 12 month in one variable
    471 !****************************************************************************************
    472           ! Get variable id
    473           CALL check_err( nf90_inq_varid(ncid, "ps", varid),"pb inq var ps" )
     531          CALL check_err(nf90_inq_varid(ncid, TRIM(cvar), varid), "pb inq var " // TRIM(cvar))
     532
    474533          ! Get the variable
    475           CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D),"pb get var ps" )
    476          
    477 ! ++) Read mass load, 12 month in one variable
    478 !****************************************************************************************
    479           ! Get variable id
    480           !CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname))
    481           IF ( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) /= nf90_noerr ) THEN
    482             WRITE(lunout,*) 'Attention load_'//TRIM(varname)//' is not in aerosol input file'
    483             load_glo2D(:,:,:)=0.0
    484           ELSE
    485             ! Get the variable
    486             CALL check_err( nf90_get_var(ncid, varid, load_glo2D),"pb get var load_"//TRIM(varname) )
    487           ENDIF
    488          
    489 ! ++) Read ap
    490 !****************************************************************************************
    491           ! Get variable id
    492           CALL check_err( nf90_inq_varid(ncid, "ap", varid),"pb inq var ap" )
    493           ! Get the variable
    494           CALL check_err( nf90_get_var(ncid, varid, pt_ap),"pb get var ap" )
    495 
    496 ! ++) Read b
    497 !****************************************************************************************
    498           ! Get variable id
    499           CALL check_err( nf90_inq_varid(ncid, "b", varid),"pb inq var b" )
    500           ! Get the variable
    501           CALL check_err( nf90_get_var(ncid, varid, pt_b),"pb get var b" )
    502 
    503          
    504 
    505        ELSE  ! old file
    506 
    507 ! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
    508 !****************************************************************************************
    509           DO imth=1, 12
    510              IF (imth==1) THEN
    511                 cvar=TRIM(varname)//'JAN'
    512              ELSE IF (imth==2) THEN
    513                 cvar=TRIM(varname)//'FEB'
    514              ELSE IF (imth==3) THEN
    515                 cvar=TRIM(varname)//'MAR'
    516              ELSE IF (imth==4) THEN
    517                 cvar=TRIM(varname)//'APR'
    518              ELSE IF (imth==5) THEN
    519                 cvar=TRIM(varname)//'MAY'
    520              ELSE IF (imth==6) THEN
    521                 cvar=TRIM(varname)//'JUN'
    522              ELSE IF (imth==7) THEN
    523                 cvar=TRIM(varname)//'JUL'
    524              ELSE IF (imth==8) THEN
    525                 cvar=TRIM(varname)//'AUG'
    526              ELSE IF (imth==9) THEN
    527                 cvar=TRIM(varname)//'SEP'
    528              ELSE IF (imth==10) THEN
    529                 cvar=TRIM(varname)//'OCT'
    530              ELSE IF (imth==11) THEN
    531                 cvar=TRIM(varname)//'NOV'
    532              ELSE IF (imth==12) THEN
    533                 cvar=TRIM(varname)//'DEC'
    534              END IF
    535              
    536              ! Get variable id
    537              CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid),"pb inq var "//TRIM(cvar) )
    538              
    539              ! Get the variable
    540              CALL check_err( nf90_get_var(ncid, varid, varmth),"pb get var "//TRIM(cvar) )
    541              
    542              ! Store in variable for the whole year
    543              varyear(:,:,:,imth)=varmth(:,:,:)
    544              
     534          CALL check_err(nf90_get_var(ncid, varid, varmth), "pb get var " // TRIM(cvar))
     535
     536          ! Store in variable for the whole year
     537          varyear(:, :, :, imth) = varmth(:, :, :)
     538
     539        END DO
     540
     541        ! Putting dummy
     542        psurf_glo2D(:, :, :) = not_valid
     543        load_glo2D(:, :, :) = not_valid
     544        pt_ap(:) = not_valid
     545        pt_b(:) = not_valid
     546
     547      END IF
     548
     549      ! 4) Close file
     550      !****************************************************************************************
     551      CALL check_err(nf90_close(ncid), "pb in close")
     552
     553
     554      ! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
     555      !****************************************************************************************
     556      ! Test if vertical levels have to be inversed
     557
     558      IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN
     559        !          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inverted'
     560        !          WRITE(lunout,*) 'before pt_ap = ', pt_ap
     561        !          WRITE(lunout,*) 'before pt_b = ', pt_b
     562
     563        ! Inverse vertical levels for varyear
     564        DO imth = 1, 12
     565          varmth(:, :, :) = varyear(:, :, :, imth) ! use varmth temporarly
     566          DO k = 1, klev_src
     567            DO j = 1, nbp_lat
     568              DO i = 1, nbp_lon
     569                varyear(i, j, k, imth) = varmth(i, j, klev_src + 1 - k)
     570              END DO
     571            END DO
    545572          END DO
    546          
    547           ! Putting dummy
    548           psurf_glo2D(:,:,:) = not_valid
    549           load_glo2D(:,:,:)  = not_valid
    550           pt_ap(:) = not_valid
    551           pt_b(:)  = not_valid
    552 
    553        END IF
    554 
    555 ! 4) Close file 
    556 !****************************************************************************************
    557        CALL check_err( nf90_close(ncid),"pb in close" )
    558      
    559 
    560 ! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
    561 !****************************************************************************************
    562 ! Test if vertical levels have to be inversed
    563 
    564        IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN
    565 !          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inverted'
    566 !          WRITE(lunout,*) 'before pt_ap = ', pt_ap
    567 !          WRITE(lunout,*) 'before pt_b = ', pt_b
    568          
    569           ! Inverse vertical levels for varyear
    570           DO imth=1, 12
    571              varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
    572              DO k=1, klev_src
    573                 DO j=1, nbp_lat
    574                    DO i=1, nbp_lon
    575                       varyear(i,j,k,imth) = varmth(i,j,klev_src+1-k)
    576                    END DO
     573        END DO
     574
     575        ! Inverte vertical axes for pt_ap and pt_b
     576        varktmp(:) = pt_ap(:)
     577        DO k = 1, klev_src
     578          pt_ap(k) = varktmp(klev_src + 1 - k)
     579        END DO
     580
     581        varktmp(:) = pt_b(:)
     582        DO k = 1, klev_src
     583          pt_b(k) = varktmp(klev_src + 1 - k)
     584        END DO
     585        WRITE(lunout, *) 'after pt_ap = ', pt_ap
     586        WRITE(lunout, *) 'after pt_b = ', pt_b
     587
     588      ELSE
     589        WRITE(lunout, *) 'Vertical axis in file ', TRIM(fname), ' is ok, no vertical inversion is done'
     590        WRITE(lunout, *) 'pt_ap = ', pt_ap
     591        WRITE(lunout, *) 'pt_b = ', pt_b
     592      END IF
     593
     594      IF (grid_type/=unstructured) THEN
     595        !     - Invert latitudes if necessary
     596        DO imth = 1, 12
     597          IF (invert_lat) THEN
     598
     599            ! Invert latitudes for the variable
     600            varmth(:, :, :) = varyear(:, :, :, imth) ! use varmth temporarly
     601            DO k = 1, klev_src
     602              DO j = 1, nbp_lat
     603                DO i = 1, nbp_lon
     604                  varyear(i, j, k, imth) = varmth(i, nbp_lat + 1 - j, k)
    577605                END DO
    578              END DO
     606              END DO
     607            END DO
     608
     609            ! Invert latitudes for surface pressure
     610            vartmp(:, :) = psurf_glo2D(:, :, imth)
     611            DO j = 1, nbp_lat
     612              DO i = 1, nbp_lon
     613                psurf_glo2D(i, j, imth) = vartmp(i, nbp_lat + 1 - j)
     614              END DO
     615            END DO
     616
     617            ! Invert latitudes for the load
     618            vartmp(:, :) = load_glo2D(:, :, imth)
     619            DO j = 1, nbp_lat
     620              DO i = 1, nbp_lon
     621                load_glo2D(i, j, imth) = vartmp(i, nbp_lat + 1 - j)
     622              END DO
     623            END DO
     624          END IF ! invert_lat
     625
     626          ! Do zonal mead at poles and distribut at whole first and last latitude
     627          DO k = 1, klev_src
     628            npole = 0.  ! North pole, j=1
     629            spole = 0.  ! South pole, j=nbp_lat
     630            DO i = 1, nbp_lon
     631              npole = npole + varyear(i, 1, k, imth)
     632              spole = spole + varyear(i, nbp_lat, k, imth)
     633            END DO
     634            npole = npole / REAL(nbp_lon)
     635            spole = spole / REAL(nbp_lon)
     636            varyear(:, 1, k, imth) = npole
     637            varyear(:, nbp_lat, k, imth) = spole
    579638          END DO
    580            
    581           ! Inverte vertical axes for pt_ap and pt_b
    582           varktmp(:) = pt_ap(:)
    583           DO k=1, klev_src
    584              pt_ap(k) = varktmp(klev_src+1-k)
    585           END DO
    586 
    587           varktmp(:) = pt_b(:)
    588           DO k=1, klev_src
    589              pt_b(k) = varktmp(klev_src+1-k)
    590           END DO
    591           WRITE(lunout,*) 'after pt_ap = ', pt_ap
    592           WRITE(lunout,*) 'after pt_b = ', pt_b
    593 
    594        ELSE
    595           WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
    596           WRITE(lunout,*) 'pt_ap = ', pt_ap
    597           WRITE(lunout,*) 'pt_b = ', pt_b
    598        END IF
    599 
    600 
    601        IF (grid_type/=unstructured) THEN
    602   !     - Invert latitudes if necessary
    603          DO imth=1, 12
    604             IF (invert_lat) THEN
    605 
    606                ! Invert latitudes for the variable
    607                varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
    608                DO k=1,klev_src
    609                   DO j=1,nbp_lat
    610                      DO i=1,nbp_lon
    611                         varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k)
    612                      END DO
    613                   END DO
    614                END DO
    615              
    616                ! Invert latitudes for surface pressure
    617                vartmp(:,:) = psurf_glo2D(:,:,imth)
    618                DO j=1,nbp_lat
    619                   DO i=1,nbp_lon
    620                      psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
    621                   END DO
    622                END DO
    623              
    624                ! Invert latitudes for the load
    625                vartmp(:,:) = load_glo2D(:,:,imth)
    626                DO j=1,nbp_lat
    627                   DO i=1,nbp_lon
    628                      load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
    629                   END DO
    630                END DO
    631             END IF ! invert_lat
    632              
    633             ! Do zonal mead at poles and distribut at whole first and last latitude
    634             DO k=1, klev_src
    635                npole=0.  ! North pole, j=1
    636                spole=0.  ! South pole, j=nbp_lat       
    637                DO i=1,nbp_lon
    638                   npole = npole + varyear(i,1,k,imth)
    639                   spole = spole + varyear(i,nbp_lat,k,imth)
    640                END DO
    641                npole = npole/REAL(nbp_lon)
    642                spole = spole/REAL(nbp_lon)
    643                varyear(:,1,    k,imth) = npole
    644                varyear(:,nbp_lat,k,imth) = spole
    645             END DO
    646          END DO ! imth
    647        
    648          ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
    649          IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1)
    650        
    651          ! Transform from 2D to 1D field
    652          CALL grid2Dto1D_glo(varyear,varyear_glo1D)
    653          CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
    654          CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
    655      
     639        END DO ! imth
     640
     641        ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat = ierr)
     642        IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3', 1)
     643
     644        ! Transform from 2D to 1D field
     645        CALL grid2Dto1D_glo(varyear, varyear_glo1D)
     646        CALL grid2Dto1D_glo(psurf_glo2D, psurf_glo1D)
     647        CALL grid2Dto1D_glo(load_glo2D, load_glo1D)
     648
    656649      ENDIF
    657      
     650
    658651    ELSE
    659         ALLOCATE(varyear_glo1D(0,0,0))       
     652      ALLOCATE(varyear_glo1D(0, 0, 0))
    660653    END IF ! is_mpi_root .AND. is_omp_root
    661654
    662 !$OMP BARRIER
    663  
    664 ! 6) Distribute to all processes
    665 !    Scatter global field(klon_glo) to local process domain(klon)
    666 !    and distribute klev_src to all processes
    667 !****************************************************************************************
     655    !$OMP BARRIER
     656
     657    ! 6) Distribute to all processes
     658    !    Scatter global field(klon_glo) to local process domain(klon)
     659    !    and distribute klev_src to all processes
     660    !****************************************************************************************
    668661
    669662    ! Distribute klev_src
    … …  
    672665    ! Allocate and distribute pt_ap and pt_b
    673666    IF (.NOT. ASSOCIATED(pt_ap)) THEN  ! if pt_ap is allocated also pt_b is allocated
    674        ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat=ierr)
    675        IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 4',1)
     667      ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat = ierr)
     668      IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 4', 1)
    676669    END IF
    677670    CALL bcast(pt_ap)
    … …  
    680673    ! Allocate space for output pointer variable at local process
    681674    IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
    682     ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
    683     ALLOCATE(pt_year_mpi(klon_mpi, klev_src, 12), stat=ierr)
    684     IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5',1)
     675    ALLOCATE(pt_year(klon, klev_src, 12), stat = ierr)
     676    ALLOCATE(pt_year_mpi(klon_mpi, klev_src, 12), stat = ierr)
     677    IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5', 1)
    685678
    686679    IF (grid_type==unstructured) THEN
    687680      IF (is_omp_master) THEN
    688         CALL xios_send_field(TRIM(varname)//"_in",varyear)
    689         CALL xios_recv_field(TRIM(varname)//"_out",pt_year_mpi)
    690         CALL xios_send_field("load_"//TRIM(varname)//"_in",load_glo2D)
    691         CALL xios_recv_field("load_"//TRIM(varname)//"_out",load_out_mpi)
     681        CALL xios_send_field(TRIM(varname) // "_in", varyear)
     682        CALL xios_recv_field(TRIM(varname) // "_out", pt_year_mpi)
     683        CALL xios_send_field("load_" // TRIM(varname) // "_in", load_glo2D)
     684        CALL xios_recv_field("load_" // TRIM(varname) // "_out", load_out_mpi)
    692685        IF (.not. allocated(psurf_interp)) THEN
    693          ! psurf_interp is a shared array
    694           ALLOCATE(psurf_interp(klon_mpi,12))
    695           CALL xios_send_field("psurf_aerosol_in",psurf_glo2D)
    696           CALL xios_recv_field("psurf_aerosol_out",psurf_interp)
     686          ! psurf_interp is a shared array
     687          ALLOCATE(psurf_interp(klon_mpi, 12))
     688          CALL xios_send_field("psurf_aerosol_in", psurf_glo2D)
     689          CALL xios_recv_field("psurf_aerosol_out", psurf_interp)
    697690        ENDIF
    698691      ENDIF
    699       CALL scatter_omp(pt_year_mpi,pt_year)
    700       CALL scatter_omp(load_out_mpi,load_out)
    701       CALL scatter_omp(psurf_interp,psurf_out)
    702       first=.FALSE.
     692      CALL scatter_omp(pt_year_mpi, pt_year)
     693      CALL scatter_omp(load_out_mpi, load_out)
     694      CALL scatter_omp(psurf_interp, psurf_out)
     695      first = .FALSE.
    703696    ELSE
    704697      ! Scatter global field to local domain at local process
    705698      CALL scatter(varyear_glo1D, pt_year)
    706699      CALL scatter(psurf_glo1D, psurf_out)
    707       CALL scatter(load_glo1D,  load_out)
     700      CALL scatter(load_glo1D, load_out)
    708701    ENDIF
    709 ! 7) Test for negative values
    710 !****************************************************************************************
     702    ! 7) Test for negative values
     703    !****************************************************************************************
    711704    IF (MINVAL(pt_year) < 0.) THEN
    712        WRITE(lunout,*) 'Warning! Negative values read from file :', fname
     705      WRITE(lunout, *) 'Warning! Negative values read from file :', fname
    713706    END IF
    714707
    … …  
    716709
    717710
    718   SUBROUTINE check_err(status,text)
     711  SUBROUTINE check_err(status, text)
    719712    USE print_control_mod, ONLY: lunout
    720713    IMPLICIT NONE
    721714
    722715    INTEGER, INTENT (IN) :: status
    723     CHARACTER(len=*), INTENT (IN), OPTIONAL :: text
     716    CHARACTER(len = *), INTENT (IN), OPTIONAL :: text
    724717
    725718    IF (status /= nf90_noerr) THEN
    726        WRITE(lunout,*) 'Error in get_aero_fromfile, netcdf error code = ',status
    727        IF (PRESENT(text)) THEN
    728           WRITE(lunout,*) 'Error in get_aero_fromfile : ',text
    729        END IF
    730        CALL abort_physic('get_aero_fromfile',trim(nf90_strerror(status)),1)
     719      WRITE(lunout, *) 'Error in get_aero_fromfile, netcdf error code = ', status
     720      IF (PRESENT(text)) THEN
     721        WRITE(lunout, *) 'Error in get_aero_fromfile : ', text
     722      END IF
     723      CALL abort_physic('get_aero_fromfile', trim(nf90_strerror(status)), 1)
    731724    END IF
    732725
Note: See TracChangeset for help on using the changeset viewer.