source: LMDZ6/branches/Amaury_dev/libf/phylmd/readaerosol_mod.F90 @ 5218

Last change on this file since 5218 was 5160, checked in by abarral, 4 months ago

Put .h into modules
File size: 30.3 KB
Line 
1! $Id: readaerosol.F90 3436 2019-01-22 16:26:21Z emillour $
2
3MODULE readaerosol_mod
4
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
11  INTEGER, SAVE :: nbp_lon_src
12  !$OMP THREADPRIVATE(nbp_lon_src)
13  INTEGER, SAVE :: nbp_lat_src
14  !$OMP THREADPRIVATE(nbp_lat_src)
15  REAL, ALLOCATABLE, SAVE :: psurf_interp(:, :)
16
17CONTAINS
18
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 lmdz_print_control, 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
94        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
95        ! pt_out has dimensions (klon, klev_src, 12)
96        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
101        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
102        ! pt_out has dimensions (klon, klev_src, 12)
103        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
104
105      ELSE
106        ! Read data from 2 decades and interpolate to actual year
107        ! a) from actual 10-yr-period
108        IF (iyr_in<1900) THEN
109          iyr1 = 1850
110          iyr2 = 1900
111        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
117        ENDIF
118
119        WRITE(cyear, '(I4)') iyr1
120        WRITE(lunout, *) 'get_aero 3 iyr_in=', iyr_in, '   ', cyear
121        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
122        ! pt_out has dimensions (klon, klev_src, 12)
123        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
124
125        ! 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
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
166        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  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")
225          ENDIF
226
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)
237          ENDIF
238          CALL xios_set_axis_attr("axis_aerosol", n_glo = klev_src)
239          CALL xios_set_fieldgroup_attr("aerosols", enabled = .TRUE.)
240
241        ENDIF
242
243      ENDIF
244    ENDIF !using_xios
245  END SUBROUTINE init_aero_fromfile
246
247
248  SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, pt_year, psurf_out, load_out)
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    !****************************************************************************************
265
266    USE dimphy
267    USE lmdz_grid_phy, ONLY: nbp_lon_ => nbp_lon, nbp_lat_ => nbp_lat, klon_glo, &
268            grid2Dto1D_glo, grid_type, unstructured
269    USE lmdz_phys_para
270    USE iophy, ONLY: io_lon, io_lat
271    USE lmdz_print_control, ONLY: lunout
272    USE lmdz_xios
273    IMPLICIT NONE
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
316    IF (grid_type==unstructured) THEN
317      nbp_lon = nbp_lon_src
318      nbp_lat = nbp_lat_src
319    ELSE
320      nbp_lon = nbp_lon_
321      nbp_lat = nbp_lat_
322    ENDIF
323
324    IF (is_mpi_root) THEN
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))
329      ALLOCATE(lon_src(nbp_lon))
330      ALLOCATE(lat_src(nbp_lat))
331      ALLOCATE(lat_src_inv(nbp_lat))
332    ELSE
333      ALLOCATE(varyear(0, 0, 0, 0))
334      ALLOCATE(psurf_glo2D(0, 0, 0))
335      ALLOCATE(load_glo2D(0, 0, 0))
336    ENDIF
337
338    ! Deallocate pointers
339    IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
340    IF (ASSOCIATED(pt_b))  DEALLOCATE(pt_b)
341
342    IF (is_mpi_root .AND. is_omp_root) THEN
343
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)
402          IF (ierr /= nf90_noerr) THEN
403            ! Dimension PRESNIVS not found either
404            CALL abort_physic('get_aero_fromfile', 'dimension lev,PRESNIVS or presnivs not in file', 1)
405          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
411          ! New file found
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
443          ierr = nf90_inq_dimid(ncid, 'time_counter', dimid)
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
530          ! Get variable id
531          CALL check_err(nf90_inq_varid(ncid, TRIM(cvar), varid), "pb inq var " // TRIM(cvar))
532
533          ! Get the variable
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
572          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)
605                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
638          END DO
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
649      ENDIF
650
651    ELSE
652      ALLOCATE(varyear_glo1D(0, 0, 0))
653    END IF ! is_mpi_root .AND. is_omp_root
654
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    !****************************************************************************************
661
662    ! Distribute klev_src
663    CALL bcast(klev_src)
664
665    ! Allocate and distribute pt_ap and pt_b
666    IF (.NOT. ASSOCIATED(pt_ap)) THEN  ! if pt_ap is allocated also pt_b is allocated
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)
669    END IF
670    CALL bcast(pt_ap)
671    CALL bcast(pt_b)
672
673    ! Allocate space for output pointer variable at local process
674    IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
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)
678
679    IF (grid_type==unstructured) THEN
680      IF (is_omp_master) THEN
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)
685        IF (.NOT. allocated(psurf_interp)) THEN
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)
690        ENDIF
691      ENDIF
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.
696    ELSE
697      ! Scatter global field to local domain at local process
698      CALL scatter(varyear_glo1D, pt_year)
699      CALL scatter(psurf_glo1D, psurf_out)
700      CALL scatter(load_glo1D, load_out)
701    ENDIF
702    ! 7) Test for negative values
703    !****************************************************************************************
704    IF (MINVAL(pt_year) < 0.) THEN
705      WRITE(lunout, *) 'Warning! Negative values read from file :', fname
706    END IF
707
708  END SUBROUTINE get_aero_fromfile
709
710
711  SUBROUTINE check_err(status, text)
712    USE lmdz_print_control, ONLY: lunout
713    IMPLICIT NONE
714
715    INTEGER, INTENT (IN) :: status
716    CHARACTER(len = *), INTENT (IN), OPTIONAL :: text
717
718    IF (status /= nf90_noerr) THEN
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)
724    END IF
725
726  END SUBROUTINE check_err
727
728
729END MODULE readaerosol_mod
Note: See TracBrowser for help on using the repository browser.