source: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/readaerosol.F90 @ 5160

Last change on this file since 5160 was 3413, checked in by Laurent Fairhead, 6 years ago

Inclusion of Yann's latest (summer/fall 2018) modifications for
convergence of DYNAMICO/LMDZ physics
YM/LF

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.3 KB
Line 
1! $Id: readaerosol.F90 3413 2018-11-12 12:52:29Z abarral $
2!
3MODULE readaerosol_mod
4
5  REAL, SAVE :: not_valid=-333.
6 
7  INTEGER, SAVE :: nbp_lon_src
8!$OMP THREADPRIVATE(nbp_lon_src) 
9  INTEGER, SAVE :: nbp_lat_src
10!$OMP THREADPRIVATE(nbp_lat_src) 
11  REAL, ALLOCATABLE, SAVE    :: psurf_interp(:,:)
12!$OMP THREADPRIVATE(psurf_interp) 
13
14CONTAINS
15
16SUBROUTINE readaerosol(name_aero, type, filename, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
17
18!****************************************************************************************
19! This routine will read the aersosol from file.
20!
21! Read a year data with get_aero_fromfile depending on aer_type :
22! - actuel   : read year 1980
23! - preind   : read natural data
24! - scenario : read one or two years and do eventually linare time interpolation
25!
26! Return pointer, pt_out, to the year read or result from interpolation
27!****************************************************************************************
28  USE dimphy
29  USE print_control_mod, ONLY: lunout
30
31  IMPLICIT NONE
32
33  ! Input arguments
34  CHARACTER(len=7), INTENT(IN) :: name_aero
35  CHARACTER(len=*), INTENT(IN) :: type  ! actuel, annuel, scenario or preind
36  CHARACTER(len=8), INTENT(IN) :: filename
37  INTEGER, INTENT(IN)          :: iyr_in
38
39  ! Output
40  INTEGER, INTENT(OUT)            :: klev_src
41  REAL, POINTER, DIMENSION(:)     :: pt_ap        ! Pointer for describing the vertical levels     
42  REAL, POINTER, DIMENSION(:)     :: pt_b         ! Pointer for describing the vertical levels     
43  REAL, POINTER, DIMENSION(:,:,:) :: pt_out       ! The massvar distributions, DIMENSION(klon, klev_src, 12)
44  REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf  ! Surface pression for 12 months
45  REAL, DIMENSION(klon,12), INTENT(OUT) :: load   ! Aerosol mass load in each column for 12 months
46
47  ! Local variables
48  CHARACTER(len=4)                :: cyear
49  REAL, POINTER, DIMENSION(:,:,:) :: pt_2
50  REAL, DIMENSION(klon,12)        :: psurf2, load2
51  INTEGER                         :: iyr1, iyr2, klev_src2
52  INTEGER                         :: it, k, i
53  LOGICAL, PARAMETER              :: lonlyone=.FALSE.
54
55!****************************************************************************************
56! Read data depending on aer_type
57!
58!****************************************************************************************
59
60  IF (type == 'actuel') THEN
61! Read and return data for year 1980
62!****************************************************************************************
63     cyear='1980'
64     ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
65     ! pt_out has dimensions (klon, klev_src, 12)
66     CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
67     
68
69  ELSE IF (type == 'preind') THEN
70! Read and return data from file with suffix .nat
71!****************************************************************************************     
72     cyear='.nat'
73     ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
74     ! pt_out has dimensions (klon, klev_src, 12)
75     CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
76     
77  ELSE IF (type == 'annuel') THEN
78! Read and return data from scenario annual files
79!****************************************************************************************     
80     WRITE(cyear,'(I4)') iyr_in
81     WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,'   ',cyear
82     ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month
83     ! pt_out has dimensions (klon, klev_src, 12)
84     CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
85     
86  ELSE IF (type == 'scenario') THEN
87! Read data depending on actual year and interpolate if necessary
88!****************************************************************************************
89     IF (iyr_in .LT. 1850) THEN
90        cyear='.nat'
91        WRITE(lunout,*) 'get_aero 1 iyr_in=', iyr_in,'   ',cyear
92        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
93        ! pt_out has dimensions (klon, klev_src, 12)
94        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
95       
96     ELSE IF (iyr_in .GE. 2100) THEN
97        cyear='2100'
98        WRITE(lunout,*) 'get_aero 2 iyr_in=', iyr_in,'   ',cyear
99        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
100        ! pt_out has dimensions (klon, klev_src, 12)
101        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
102       
103     ELSE
104        ! Read data from 2 decades and interpolate to actual year
105        ! a) from actual 10-yr-period
106        IF (iyr_in.LT.1900) THEN
107           iyr1 = 1850
108           iyr2 = 1900
109        ELSE IF (iyr_in.GE.1900.AND.iyr_in.LT.1920) THEN
110           iyr1 = 1900
111           iyr2 = 1920
112        ELSE
113           iyr1 = INT(iyr_in/10)*10
114           iyr2 = INT(1+iyr_in/10)*10
115        ENDIF
116       
117        WRITE(cyear,'(I4)') iyr1
118        WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,'   ',cyear
119        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
120        ! pt_out has dimensions (klon, klev_src, 12)
121        CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
122       
123        ! If to read two decades:
124        IF (.NOT.lonlyone) THEN
125           
126           ! b) from the next following one
127           WRITE(cyear,'(I4)') iyr2
128           WRITE(lunout,*) 'get_aero 4 iyr_in=', iyr_in,'   ',cyear
129           
130           NULLIFY(pt_2)
131           ! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month
132           ! pt_2 has dimensions (klon, klev_src, 12)
133           CALL get_aero_fromfile(name_aero, cyear, filename, klev_src2, pt_ap, pt_b, pt_2, psurf2, load2)
134           ! Test for same number of vertical levels
135           IF (klev_src /= klev_src2) THEN
136              WRITE(lunout,*) 'Two aerosols files with different number of vertical levels is not allowded'
137              CALL abort_physic('readaersosol','Error in number of vertical levels',1)
138           END IF
139           
140           ! Linare interpolate to the actual year:
141           DO it=1,12
142              DO k=1,klev_src
143                 DO i = 1, klon
144                    pt_out(i,k,it) = &
145                         pt_out(i,k,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * &
146                         (pt_out(i,k,it) - pt_2(i,k,it))
147                 END DO
148              END DO
149
150              DO i = 1, klon
151                 psurf(i,it) = &
152                      psurf(i,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * &
153                      (psurf(i,it) - psurf2(i,it))
154
155                 load(i,it) = &
156                      load(i,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * &
157                      (load(i,it) - load2(i,it))
158              END DO
159           END DO
160
161           ! Deallocate pt_2 no more needed
162           DEALLOCATE(pt_2)
163           
164        END IF ! lonlyone
165     END IF ! iyr_in .LT. 1850
166
167  ELSE
168     WRITE(lunout,*)'This option is not implemented : aer_type = ', type,' name_aero=',name_aero
169     CALL abort_physic('readaerosol','Error : aer_type parameter not accepted',1)
170  END IF ! type
171
172
173END SUBROUTINE readaerosol
174
175
176  SUBROUTINE init_aero_fromfile(flag_aerosol)
177  USE netcdf
178  USE mod_phys_lmdz_para
179  USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured
180  USE xios
181  IMPLICIT NONE
182  INTEGER, INTENT(IN) :: flag_aerosol
183  REAL,ALLOCATABLE :: lat_src(:)
184  REAL,ALLOCATABLE :: lon_src(:)
185  CHARACTER(LEN=*),PARAMETER :: file_aerosol='aerosols.nat.nc'
186  CHARACTER(LEN=*),PARAMETER :: file_so4='so4.nat.nc'
187  INTEGER :: klev_src
188  INTEGER :: ierr ,ncid, dimID, varid
189  REAL :: null_array(0)
190
191  IF (flag_aerosol>0 .AND. grid_type==unstructured) THEN
192 
193    IF (is_omp_root) THEN
194 
195      IF (is_mpi_root) THEN
196   
197        IF (nf90_open(TRIM(file_aerosol), NF90_NOWRITE, ncid) /= NF90_NOERR) THEN
198          CALL check_err( nf90_open(TRIM(file_so4), NF90_NOWRITE, ncid), "pb open "//trim(file_so4) )
199        ENDIF
200
201        ! Read and test longitudes
202        CALL check_err( nf90_inq_dimid(ncid, "lon", dimID),"pb inq dim lon")
203        CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lon_src),"pb inq dim lon")
204        CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
205        ALLOCATE(lon_src(nbp_lon_src))
206        CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
207
208        ! Read and test latitudes
209        CALL check_err( nf90_inq_dimid(ncid, "lat", dimID),"pb inq dim lat")
210        CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lat_src),"pb inq dim lat")
211        CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
212        ALLOCATE(lat_src(nbp_lat_src))
213        CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
214        IF (nf90_inq_dimid(ncid, 'lev', dimid) /= NF90_NOERR) THEN
215          IF (nf90_inq_dimid(ncid, 'presnivs', dimid)/= NF90_NOERR) THEN
216             CALL check_err(nf90_inq_dimid(ncid, 'PRESNIVS', dimid),'dimension lev,PRESNIVS or presnivs not in file')
217          ENDIF
218        ENDIF
219        CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
220        CALL check_err( nf90_close(ncid),"pb in close" )   
221      ENDIF
222
223      CALL bcast_mpi(nbp_lat_src)
224      CALL bcast_mpi(nbp_lon_src)
225      CALL bcast_mpi(klev_src)
226
227      IF (is_mpi_root ) THEN
228        CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=nbp_lat_src, jbegin=0, latvalue_1d=lat_src)
229        CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=nbp_lon_src, ibegin=0, lonvalue_1d=lon_src)
230      ELSE
231        CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=0, jbegin=0, latvalue_1d=null_array )
232        CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=0, ibegin=0, lonvalue_1d=null_array)
233      ENDIF
234      CALL xios_set_axis_attr("axis_aerosol",n_glo=klev_src)
235      CALL xios_set_fieldgroup_attr("aerosols", enabled=.TRUE.)
236 
237    ENDIF
238   
239  ENDIF   
240 
241  END SUBROUTINE init_aero_fromfile
242
243
244   
245
246  SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, pt_year, psurf_out, load_out)
247!****************************************************************************************
248! Read 12 month aerosol from file and distribute to local process on physical grid.
249! Vertical levels, klev_src, may differ from model levels if new file format.
250!
251! For mpi_root and master thread :
252! 1) Open file
253! 2) Find vertical dimension klev_src
254! 3) Read field month by month
255! 4) Close file 
256! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
257!     - Also the levels and the latitudes have to be inversed
258!
259! For all processes and threads :
260! 6) Scatter global field(klon_glo) to local process domain(klon)
261! 7) Test for negative values
262!****************************************************************************************
263
264    USE netcdf
265    USE dimphy
266    USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, klon_glo, &
267                                 grid2Dto1D_glo, grid_type, unstructured
268    USE mod_phys_lmdz_para
269    USE iophy, ONLY : io_lon, io_lat
270    USE print_control_mod, ONLY: lunout
271    USE xios
272
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
353       IF (grid_type/=unstructured) THEN
354
355! Test for equal longitudes and latitudes in file and model
356!****************************************************************************************
357       ! Read and test longitudes
358         CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
359         CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
360       
361         IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
362            WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
363            WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
364            WRITE(lunout,*) 'longitudes in model :', io_lon
365         
366            CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
367         END IF
368
369         ! Read and test latitudes
370         CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
371         CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
372
373         ! Invert source latitudes
374         DO j = 1, nbp_lat
375            lat_src_inv(j) = lat_src(nbp_lat +1 -j)
376         END DO
377
378         IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
379            ! Latitudes are the same
380            invert_lat=.FALSE.
381         ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
382            ! Inverted source latitudes correspond to model latitudes
383            WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
384            invert_lat=.TRUE.
385         ELSE
386            WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
387            WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
388            WRITE(lunout,*) 'latitudes in model :', io_lat
389            CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
390         END IF
391       ENDIF
392
393! 2) Check if old or new file is avalabale.
394!    New type of file should contain the dimension 'lev'
395!    Old type of file should contain the dimension 'PRESNIVS'
396!****************************************************************************************
397       ierr = nf90_inq_dimid(ncid, 'lev', dimid)
398       IF (ierr /= NF90_NOERR) THEN
399          ! Coordinate axe lev not found. Check for presnivs.
400          ierr = nf90_inq_dimid(ncid, 'presnivs', dimid)
401          IF (ierr /= NF90_NOERR) THEN
402             ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid)
403             IF (ierr /= NF90_NOERR) THEN
404                ! Dimension PRESNIVS not found either
405                CALL abort_physic('get_aero_fromfile', 'dimension lev,PRESNIVS or presnivs not in file',1)
406             ELSE
407                ! Old file found
408                new_file=.FALSE.
409                WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done'
410             END IF
411          ELSE
412             ! New file found
413             new_file=.TRUE.
414             WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
415          ENDIF
416       ELSE
417          ! New file found
418          new_file=.TRUE.
419          WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
420       END IF
421       
422! 2) Find vertical dimension klev_src
423!****************************************************************************************
424       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
425       
426     ! Allocate variables depending on the number of vertical levels
427       ALLOCATE(varmth(nbp_lon,nbp_lat, klev_src), varyear(nbp_lon,nbp_lat, klev_src, 12), stat=ierr)
428       IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 1',1)
429
430       ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr)
431       IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 2',1)
432
433! 3) Read all variables from file
434!    There is 2 options for the file structure :
435!    new_file=TRUE  : read varyear, ps, pt_ap and pt_b
436!    new_file=FALSE : read varyear month by month
437!****************************************************************************************
438
439       IF (new_file) THEN
440! ++) Check number of month in file opened
441!**************************************************************************************************
442       ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
443       if (ierr /= NF90_NOERR) THEN
444          ierr = nf90_inq_dimid(ncid, 'time_counter', dimid)
445       ENDIF
446       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps),"pb inq dim TIME or time_counter" )
447!       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
448       IF (nbr_tsteps /= 12 ) THEN
449    CALL abort_physic('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)' &
450                     ,1)
451       ENDIF
452
453! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
454!****************************************************************************************
455          ! Get variable id
456          !CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) )
457          print *,'readaerosol ', TRIM(varname)
458          IF ( nf90_inq_varid(ncid, TRIM(varname), varid) /= NF90_NOERR ) THEN
459            ! Variable is not there
460            WRITE(lunout,*) 'Attention '//TRIM(varname)//' is not in aerosol input file'
461            varyear(:,:,:,:)=0.0
462          ELSE
463            ! Get the variable
464            CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)),"pb get var "//TRIM(varname) )
465          ENDIF
466         
467! ++) Read surface pression, 12 month in one variable
468!****************************************************************************************
469          ! Get variable id
470          CALL check_err( nf90_inq_varid(ncid, "ps", varid),"pb inq var ps" )
471          ! Get the variable
472          CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D),"pb get var ps" )
473         
474! ++) Read mass load, 12 month in one variable
475!****************************************************************************************
476          ! Get variable id
477          !CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname))
478          IF ( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) /= NF90_NOERR ) THEN
479            WRITE(lunout,*) 'Attention load_'//TRIM(varname)//' is not in aerosol input file'
480            load_glo2D(:,:,:)=0.0
481          ELSE
482            ! Get the variable
483            CALL check_err( nf90_get_var(ncid, varid, load_glo2D),"pb get var load_"//TRIM(varname) )
484          ENDIF
485         
486! ++) Read ap
487!****************************************************************************************
488          ! Get variable id
489          CALL check_err( nf90_inq_varid(ncid, "ap", varid),"pb inq var ap" )
490          ! Get the variable
491          CALL check_err( nf90_get_var(ncid, varid, pt_ap),"pb get var ap" )
492
493! ++) Read b
494!****************************************************************************************
495          ! Get variable id
496          CALL check_err( nf90_inq_varid(ncid, "b", varid),"pb inq var b" )
497          ! Get the variable
498          CALL check_err( nf90_get_var(ncid, varid, pt_b),"pb get var b" )
499
500         
501
502       ELSE  ! old file
503
504! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
505!****************************************************************************************
506          DO imth=1, 12
507             IF (imth.EQ.1) THEN
508                cvar=TRIM(varname)//'JAN'
509             ELSE IF (imth.EQ.2) THEN
510                cvar=TRIM(varname)//'FEB'
511             ELSE IF (imth.EQ.3) THEN
512                cvar=TRIM(varname)//'MAR'
513             ELSE IF (imth.EQ.4) THEN
514                cvar=TRIM(varname)//'APR'
515             ELSE IF (imth.EQ.5) THEN
516                cvar=TRIM(varname)//'MAY'
517             ELSE IF (imth.EQ.6) THEN
518                cvar=TRIM(varname)//'JUN'
519             ELSE IF (imth.EQ.7) THEN
520                cvar=TRIM(varname)//'JUL'
521             ELSE IF (imth.EQ.8) THEN
522                cvar=TRIM(varname)//'AUG'
523             ELSE IF (imth.EQ.9) THEN
524                cvar=TRIM(varname)//'SEP'
525             ELSE IF (imth.EQ.10) THEN
526                cvar=TRIM(varname)//'OCT'
527             ELSE IF (imth.EQ.11) THEN
528                cvar=TRIM(varname)//'NOV'
529             ELSE IF (imth.EQ.12) THEN
530                cvar=TRIM(varname)//'DEC'
531             END IF
532             
533             ! Get variable id
534             CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid),"pb inq var "//TRIM(cvar) )
535             
536             ! Get the variable
537             CALL check_err( nf90_get_var(ncid, varid, varmth),"pb get var "//TRIM(cvar) )
538             
539             ! Store in variable for the whole year
540             varyear(:,:,:,imth)=varmth(:,:,:)
541             
542          END DO
543         
544          ! Putting dummy
545          psurf_glo2D(:,:,:) = not_valid
546          load_glo2D(:,:,:)  = not_valid
547          pt_ap(:) = not_valid
548          pt_b(:)  = not_valid
549
550       END IF
551
552! 4) Close file 
553!****************************************************************************************
554       CALL check_err( nf90_close(ncid),"pb in close" )
555     
556
557! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
558!****************************************************************************************
559! Test if vertical levels have to be inversed
560
561       IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN
562!          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inverted'
563!          WRITE(lunout,*) 'before pt_ap = ', pt_ap
564!          WRITE(lunout,*) 'before pt_b = ', pt_b
565         
566          ! Inverse vertical levels for varyear
567          DO imth=1, 12
568             varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
569             DO k=1, klev_src
570                DO j=1, nbp_lat
571                   DO i=1, nbp_lon
572                      varyear(i,j,k,imth) = varmth(i,j,klev_src+1-k)
573                   END DO
574                END DO
575             END DO
576          END DO
577           
578          ! Inverte vertical axes for pt_ap and pt_b
579          varktmp(:) = pt_ap(:)
580          DO k=1, klev_src
581             pt_ap(k) = varktmp(klev_src+1-k)
582          END DO
583
584          varktmp(:) = pt_b(:)
585          DO k=1, klev_src
586             pt_b(k) = varktmp(klev_src+1-k)
587          END DO
588          WRITE(lunout,*) 'after pt_ap = ', pt_ap
589          WRITE(lunout,*) 'after pt_b = ', pt_b
590
591       ELSE
592          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
593          WRITE(lunout,*) 'pt_ap = ', pt_ap
594          WRITE(lunout,*) 'pt_b = ', pt_b
595       END IF
596
597
598       IF (grid_type/=unstructured) THEN
599  !     - Invert latitudes if necessary
600         DO imth=1, 12
601            IF (invert_lat) THEN
602
603               ! Invert latitudes for the variable
604               varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
605               DO k=1,klev_src
606                  DO j=1,nbp_lat
607                     DO i=1,nbp_lon
608                        varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k)
609                     END DO
610                  END DO
611               END DO
612             
613               ! Invert latitudes for surface pressure
614               vartmp(:,:) = psurf_glo2D(:,:,imth)
615               DO j=1,nbp_lat
616                  DO i=1,nbp_lon
617                     psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
618                  END DO
619               END DO
620             
621               ! Invert latitudes for the load
622               vartmp(:,:) = load_glo2D(:,:,imth)
623               DO j=1,nbp_lat
624                  DO i=1,nbp_lon
625                     load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
626                  END DO
627               END DO
628            END IF ! invert_lat
629             
630            ! Do zonal mead at poles and distribut at whole first and last latitude
631            DO k=1, klev_src
632               npole=0.  ! North pole, j=1
633               spole=0.  ! South pole, j=nbp_lat       
634               DO i=1,nbp_lon
635                  npole = npole + varyear(i,1,k,imth)
636                  spole = spole + varyear(i,nbp_lat,k,imth)
637               END DO
638               npole = npole/REAL(nbp_lon)
639               spole = spole/REAL(nbp_lon)
640               varyear(:,1,    k,imth) = npole
641               varyear(:,nbp_lat,k,imth) = spole
642            END DO
643         END DO ! imth
644       
645         ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
646         IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1)
647       
648         ! Transform from 2D to 1D field
649         CALL grid2Dto1D_glo(varyear,varyear_glo1D)
650         CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
651         CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
652     
653      ENDIF
654     
655    ELSE
656        ALLOCATE(varyear_glo1D(0,0,0))       
657    END IF ! is_mpi_root .AND. is_omp_root
658
659!$OMP BARRIER
660 
661! 6) Distribute to all processes
662!    Scatter global field(klon_glo) to local process domain(klon)
663!    and distribute klev_src to all processes
664!****************************************************************************************
665
666    ! Distribute klev_src
667    CALL bcast(klev_src)
668
669    ! Allocate and distribute pt_ap and pt_b
670    IF (.NOT. ASSOCIATED(pt_ap)) THEN  ! if pt_ap is allocated also pt_b is allocated
671       ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat=ierr)
672       IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 4',1)
673    END IF
674    CALL bcast(pt_ap)
675    CALL bcast(pt_b)
676
677    ! Allocate space for output pointer variable at local process
678    IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
679    ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
680    ALLOCATE(pt_year_mpi(klon_mpi, klev_src, 12), stat=ierr)
681    IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5',1)
682
683    IF (grid_type==unstructured) THEN
684      IF (is_omp_master) THEN
685        CALL xios_send_field(TRIM(varname)//"_in",varyear)
686        CALL xios_recv_field(TRIM(varname)//"_out",pt_year_mpi)
687        CALL xios_send_field("load_"//TRIM(varname)//"_in",load_glo2D)
688        CALL xios_recv_field("load_"//TRIM(varname)//"_out",load_out_mpi)
689        IF (first) THEN
690          ALLOCATE(psurf_interp(klon_mpi,12))
691          CALL xios_send_field("psurf_aerosol_in",psurf_glo2D)
692          CALL xios_recv_field("psurf_aerosol_out",psurf_interp)
693        ENDIF
694      ENDIF
695      CALL scatter_omp(pt_year_mpi,pt_year)
696      CALL scatter_omp(load_out_mpi,load_out)
697      CALL scatter_omp(psurf_interp,psurf_out)
698      first=.FALSE.
699    ELSE
700      ! Scatter global field to local domain at local process
701      CALL scatter(varyear_glo1D, pt_year)
702      CALL scatter(psurf_glo1D, psurf_out)
703      CALL scatter(load_glo1D,  load_out)
704    ENDIF
705! 7) Test for negative values
706!****************************************************************************************
707    IF (MINVAL(pt_year) < 0.) THEN
708       WRITE(lunout,*) 'Warning! Negative values read from file :', fname
709    END IF
710
711  END SUBROUTINE get_aero_fromfile
712
713
714  SUBROUTINE check_err(status,text)
715    USE netcdf
716    USE print_control_mod, ONLY: lunout
717    IMPLICIT NONE
718
719    INTEGER, INTENT (IN) :: status
720    CHARACTER(len=*), INTENT (IN), OPTIONAL :: text
721
722    IF (status /= NF90_NOERR) THEN
723       WRITE(lunout,*) 'Error in get_aero_fromfile, netcdf error code = ',status
724       IF (PRESENT(text)) THEN
725          WRITE(lunout,*) 'Error in get_aero_fromfile : ',text
726       END IF
727       CALL abort_physic('get_aero_fromfile',trim(nf90_strerror(status)),1)
728    END IF
729
730  END SUBROUTINE check_err
731
732
733END MODULE readaerosol_mod
Note: See TracBrowser for help on using the repository browser.