source: LMDZ6/trunk/libf/phylmd/readaerosol_mod.F90 @ 3847

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

For checking purposes

File size: 30.4 KB
RevLine 
[3440]1! $Id: readaerosol.F90 3436 2019-01-22 16:26:21Z emillour $
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#ifdef CPP_XIOS
181  USE xios
182#endif
183  IMPLICIT NONE
184  INTEGER, INTENT(IN) :: flag_aerosol
185#ifdef CPP_XIOS
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 (flag_aerosol>0 .AND. grid_type==unstructured) THEN
195 
196    IF (is_omp_root) THEN
197 
198      IF (is_mpi_root) THEN
199   
200        IF (nf90_open(TRIM(file_aerosol), NF90_NOWRITE, ncid) /= NF90_NOERR) THEN
201          CALL check_err( nf90_open(TRIM(file_so4), NF90_NOWRITE, ncid), "pb open "//trim(file_so4) )
202        ENDIF
203
204        ! Read and test longitudes
205        CALL check_err( nf90_inq_dimid(ncid, "lon", dimID),"pb inq dim lon")
206        CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lon_src),"pb inq dim lon")
207        CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
208        ALLOCATE(lon_src(nbp_lon_src))
209        CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
210
211        ! Read and test latitudes
212        CALL check_err( nf90_inq_dimid(ncid, "lat", dimID),"pb inq dim lat")
213        CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lat_src),"pb inq dim lat")
214        CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
215        ALLOCATE(lat_src(nbp_lat_src))
216        CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
217        IF (nf90_inq_dimid(ncid, 'lev', dimid) /= NF90_NOERR) THEN
218          IF (nf90_inq_dimid(ncid, 'presnivs', dimid)/= NF90_NOERR) THEN
219             CALL check_err(nf90_inq_dimid(ncid, 'PRESNIVS', dimid),'dimension lev,PRESNIVS or presnivs not in file')
220          ENDIF
221        ENDIF
222        CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
223        CALL check_err( nf90_close(ncid),"pb in close" )   
224      ENDIF
225
226      CALL bcast_mpi(nbp_lat_src)
227      CALL bcast_mpi(nbp_lon_src)
228      CALL bcast_mpi(klev_src)
229
230      IF (is_mpi_root ) THEN
231        CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=nbp_lat_src, jbegin=0, latvalue_1d=lat_src)
232        CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=nbp_lon_src, ibegin=0, lonvalue_1d=lon_src)
233      ELSE
234        CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=0, jbegin=0, latvalue_1d=null_array )
235        CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=0, ibegin=0, lonvalue_1d=null_array)
236      ENDIF
237      CALL xios_set_axis_attr("axis_aerosol",n_glo=klev_src)
238      CALL xios_set_fieldgroup_attr("aerosols", enabled=.TRUE.)
239 
240    ENDIF
241   
242  ENDIF   
243#endif
244  END SUBROUTINE init_aero_fromfile
245
246
247   
248
249  SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, pt_year, psurf_out, load_out)
250!****************************************************************************************
251! Read 12 month aerosol from file and distribute to local process on physical grid.
252! Vertical levels, klev_src, may differ from model levels if new file format.
253!
254! For mpi_root and master thread :
255! 1) Open file
256! 2) Find vertical dimension klev_src
257! 3) Read field month by month
258! 4) Close file 
259! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
260!     - Also the levels and the latitudes have to be inversed
261!
262! For all processes and threads :
263! 6) Scatter global field(klon_glo) to local process domain(klon)
264! 7) Test for negative values
265!****************************************************************************************
266
267    USE netcdf
268    USE dimphy
269    USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, klon_glo, &
270                                 grid2Dto1D_glo, grid_type, unstructured
271    USE mod_phys_lmdz_para
272    USE iophy, ONLY : io_lon, io_lat
273    USE print_control_mod, ONLY: lunout
274#ifdef CPP_XIOS
275    USE xios
276#endif
277    IMPLICIT NONE
278     
279! Input argumets
280    CHARACTER(len=7), INTENT(IN)          :: varname
281    CHARACTER(len=4), INTENT(IN)          :: cyr
282    CHARACTER(len=8), INTENT(IN)          :: filename
283
284! Output arguments
285    INTEGER, INTENT(OUT)                  :: klev_src     ! Number of vertical levels in file
286    REAL, POINTER, DIMENSION(:)           :: pt_ap        ! Pointer for describing the vertical levels     
287    REAL, POINTER, DIMENSION(:)           :: pt_b         ! Pointer for describing the vertical levels     
288    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year      ! Pointer-variabale from file, 12 month, grid : klon,klev_src
289    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year_mpi  ! Pointer-variabale from file, 12 month, grid : klon,klev_src
290    REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out    ! Surface pression for 12 months
291    REAL, DIMENSION(klon_mpi,12)          :: psurf_out_mpi    ! Surface pression for 12 months
292    REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out     ! Aerosol mass load in each column
293    REAL, DIMENSION(klon_mpi,12)          :: load_out_mpi     ! Aerosol mass load in each column
294    INTEGER                               :: nbr_tsteps   ! number of month in file read
295
296! Local variables
297    CHARACTER(len=30)     :: fname
298    CHARACTER(len=30)     :: cvar
299    INTEGER               :: ncid, dimid, varid
300    INTEGER               :: imth, i, j, k, ierr
301    REAL                  :: npole, spole
302    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varmth
303    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear       ! Global variable read from file, 12 month
304    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varyear_glo1D !(klon_glo, klev_src, 12)
305    REAL, ALLOCATABLE, DIMENSION(:)       :: varktmp
306
307    REAL,  ALLOCATABLE                    :: psurf_glo2D(:,:,:)   ! Surface pression for 12 months on dynamics global grid
308    REAL, DIMENSION(klon_glo,12)          :: psurf_glo1D   ! -"- on physical global grid
309    REAL,  ALLOCATABLE                    :: load_glo2D(:,:,:)    ! Load for 12 months on dynamics global grid
310    REAL, DIMENSION(klon_glo,12)          :: load_glo1D    ! -"- on physical global grid
311    REAL, ALLOCATABLE, DIMENSION(:,:)     :: vartmp
312    REAL, ALLOCATABLE,DIMENSION(:)        :: lon_src              ! longitudes in file
313    REAL, ALLOCATABLE,DIMENSION(:)        :: lat_src, lat_src_inv ! latitudes in file
314    LOGICAL                               :: new_file             ! true if new file format detected
315    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
316    INTEGER                               :: nbp_lon, nbp_lat
317    LOGICAL,SAVE                          :: first=.TRUE.
318!$OMP THREADPRIVATE(first)
319   
320    IF (grid_type==unstructured) THEN
321      nbp_lon=nbp_lon_src
322      nbp_lat=nbp_lat_src
323    ELSE
324      nbp_lon=nbp_lon_
325      nbp_lat=nbp_lat_
326    ENDIF
327   
328    IF (is_mpi_root) THEN
329   
330      ALLOCATE(psurf_glo2D(nbp_lon,nbp_lat,12))
331      ALLOCATE(load_glo2D(nbp_lon,nbp_lat,12))
332      ALLOCATE(vartmp(nbp_lon,nbp_lat))
333      ALLOCATE(lon_src(nbp_lon))
334      ALLOCATE(lat_src(nbp_lat))
335      ALLOCATE(lat_src_inv(nbp_lat))
336    ELSE
337      ALLOCATE(varyear(0,0,0,0))
338      ALLOCATE(psurf_glo2D(0,0,0))
339      ALLOCATE(load_glo2D(0,0,0))
340    ENDIF
341           
342    ! Deallocate pointers
343    IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
344    IF (ASSOCIATED(pt_b))  DEALLOCATE(pt_b)
345
346    IF (is_mpi_root .AND. is_omp_root) THEN
347
348! 1) Open file
349!****************************************************************************************
350! Add suffix to filename
351       fname = trim(filename)//cyr//'.nc'
352 
353       WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname)
354       CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(fname) )
355
356
357       IF (grid_type/=unstructured) THEN
358
359! Test for equal longitudes and latitudes in file and model
360!****************************************************************************************
361       ! Read and test longitudes
362         CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
363         CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
364       
365         IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
366            WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
367            WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
368            WRITE(lunout,*) 'longitudes in model :', io_lon
369         
370            CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
371         END IF
372
373         ! Read and test latitudes
374         CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
375         CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
376
377         ! Invert source latitudes
378         DO j = 1, nbp_lat
379            lat_src_inv(j) = lat_src(nbp_lat +1 -j)
380         END DO
381
382         IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
383            ! Latitudes are the same
384            invert_lat=.FALSE.
385         ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
386            ! Inverted source latitudes correspond to model latitudes
387            WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
388            invert_lat=.TRUE.
389         ELSE
390            WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
391            WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
392            WRITE(lunout,*) 'latitudes in model :', io_lat
393            CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
394         END IF
395       ENDIF
396
397! 2) Check if old or new file is avalabale.
398!    New type of file should contain the dimension 'lev'
399!    Old type of file should contain the dimension 'PRESNIVS'
400!****************************************************************************************
401       ierr = nf90_inq_dimid(ncid, 'lev', dimid)
402       IF (ierr /= NF90_NOERR) THEN
403          ! Coordinate axe lev not found. Check for presnivs.
404          ierr = nf90_inq_dimid(ncid, 'presnivs', dimid)
405          IF (ierr /= NF90_NOERR) THEN
406             ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid)
407             IF (ierr /= NF90_NOERR) THEN
408                ! Dimension PRESNIVS not found either
409                CALL abort_physic('get_aero_fromfile', 'dimension lev,PRESNIVS or presnivs not in file',1)
410             ELSE
411                ! Old file found
412                new_file=.FALSE.
413                WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done'
414             END IF
415          ELSE
416             ! New file found
417             new_file=.TRUE.
418             WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
419          ENDIF
420       ELSE
421          ! New file found
422          new_file=.TRUE.
423          WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
424       END IF
425       
426! 2) Find vertical dimension klev_src
427!****************************************************************************************
428       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
429       
430     ! Allocate variables depending on the number of vertical levels
431       ALLOCATE(varmth(nbp_lon,nbp_lat, klev_src), varyear(nbp_lon,nbp_lat, klev_src, 12), stat=ierr)
432       IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 1',1)
433
434       ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr)
435       IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 2',1)
436
437! 3) Read all variables from file
438!    There is 2 options for the file structure :
439!    new_file=TRUE  : read varyear, ps, pt_ap and pt_b
440!    new_file=FALSE : read varyear month by month
441!****************************************************************************************
442
443       IF (new_file) THEN
444! ++) Check number of month in file opened
445!**************************************************************************************************
446       ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
447       if (ierr /= NF90_NOERR) THEN
448          ierr = nf90_inq_dimid(ncid, 'time_counter', dimid)
449       ENDIF
450       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps),"pb inq dim TIME or time_counter" )
451!       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
452       IF (nbr_tsteps /= 12 ) THEN
453    CALL abort_physic('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)' &
454                     ,1)
455       ENDIF
456
457! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
458!****************************************************************************************
459          ! Get variable id
460          !CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) )
461          print *,'readaerosol ', TRIM(varname)
462          IF ( nf90_inq_varid(ncid, TRIM(varname), varid) /= NF90_NOERR ) THEN
463            ! Variable is not there
464            WRITE(lunout,*) 'Attention '//TRIM(varname)//' is not in aerosol input file'
465            varyear(:,:,:,:)=0.0
466          ELSE
467            ! Get the variable
468            CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)),"pb get var "//TRIM(varname) )
469          ENDIF
470         
471! ++) Read surface pression, 12 month in one variable
472!****************************************************************************************
473          ! Get variable id
474          CALL check_err( nf90_inq_varid(ncid, "ps", varid),"pb inq var ps" )
475          ! Get the variable
476          CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D),"pb get var ps" )
477         
478! ++) Read mass load, 12 month in one variable
479!****************************************************************************************
480          ! Get variable id
481          !CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname))
482          IF ( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) /= NF90_NOERR ) THEN
483            WRITE(lunout,*) 'Attention load_'//TRIM(varname)//' is not in aerosol input file'
484            load_glo2D(:,:,:)=0.0
485          ELSE
486            ! Get the variable
487            CALL check_err( nf90_get_var(ncid, varid, load_glo2D),"pb get var load_"//TRIM(varname) )
488          ENDIF
489         
490! ++) Read ap
491!****************************************************************************************
492          ! Get variable id
493          CALL check_err( nf90_inq_varid(ncid, "ap", varid),"pb inq var ap" )
494          ! Get the variable
495          CALL check_err( nf90_get_var(ncid, varid, pt_ap),"pb get var ap" )
496
497! ++) Read b
498!****************************************************************************************
499          ! Get variable id
500          CALL check_err( nf90_inq_varid(ncid, "b", varid),"pb inq var b" )
501          ! Get the variable
502          CALL check_err( nf90_get_var(ncid, varid, pt_b),"pb get var b" )
503
504         
505
506       ELSE  ! old file
507
508! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
509!****************************************************************************************
510          DO imth=1, 12
511             IF (imth.EQ.1) THEN
512                cvar=TRIM(varname)//'JAN'
513             ELSE IF (imth.EQ.2) THEN
514                cvar=TRIM(varname)//'FEB'
515             ELSE IF (imth.EQ.3) THEN
516                cvar=TRIM(varname)//'MAR'
517             ELSE IF (imth.EQ.4) THEN
518                cvar=TRIM(varname)//'APR'
519             ELSE IF (imth.EQ.5) THEN
520                cvar=TRIM(varname)//'MAY'
521             ELSE IF (imth.EQ.6) THEN
522                cvar=TRIM(varname)//'JUN'
523             ELSE IF (imth.EQ.7) THEN
524                cvar=TRIM(varname)//'JUL'
525             ELSE IF (imth.EQ.8) THEN
526                cvar=TRIM(varname)//'AUG'
527             ELSE IF (imth.EQ.9) THEN
528                cvar=TRIM(varname)//'SEP'
529             ELSE IF (imth.EQ.10) THEN
530                cvar=TRIM(varname)//'OCT'
531             ELSE IF (imth.EQ.11) THEN
532                cvar=TRIM(varname)//'NOV'
533             ELSE IF (imth.EQ.12) THEN
534                cvar=TRIM(varname)//'DEC'
535             END IF
536             
537             ! Get variable id
538             CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid),"pb inq var "//TRIM(cvar) )
539             
540             ! Get the variable
541             CALL check_err( nf90_get_var(ncid, varid, varmth),"pb get var "//TRIM(cvar) )
542             
543             ! Store in variable for the whole year
544             varyear(:,:,:,imth)=varmth(:,:,:)
545             
546          END DO
547         
548          ! Putting dummy
549          psurf_glo2D(:,:,:) = not_valid
550          load_glo2D(:,:,:)  = not_valid
551          pt_ap(:) = not_valid
552          pt_b(:)  = not_valid
553
554       END IF
555
556! 4) Close file 
557!****************************************************************************************
558       CALL check_err( nf90_close(ncid),"pb in close" )
559     
560
561! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
562!****************************************************************************************
563! Test if vertical levels have to be inversed
564
565       IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN
566!          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inverted'
567!          WRITE(lunout,*) 'before pt_ap = ', pt_ap
568!          WRITE(lunout,*) 'before pt_b = ', pt_b
569         
570          ! Inverse vertical levels for varyear
571          DO imth=1, 12
572             varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
573             DO k=1, klev_src
574                DO j=1, nbp_lat
575                   DO i=1, nbp_lon
576                      varyear(i,j,k,imth) = varmth(i,j,klev_src+1-k)
577                   END DO
578                END DO
579             END DO
580          END DO
581           
582          ! Inverte vertical axes for pt_ap and pt_b
583          varktmp(:) = pt_ap(:)
584          DO k=1, klev_src
585             pt_ap(k) = varktmp(klev_src+1-k)
586          END DO
587
588          varktmp(:) = pt_b(:)
589          DO k=1, klev_src
590             pt_b(k) = varktmp(klev_src+1-k)
591          END DO
592          WRITE(lunout,*) 'after pt_ap = ', pt_ap
593          WRITE(lunout,*) 'after pt_b = ', pt_b
594
595       ELSE
596          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
597          WRITE(lunout,*) 'pt_ap = ', pt_ap
598          WRITE(lunout,*) 'pt_b = ', pt_b
599       END IF
600
601
602       IF (grid_type/=unstructured) THEN
603  !     - Invert latitudes if necessary
604         DO imth=1, 12
605            IF (invert_lat) THEN
606
607               ! Invert latitudes for the variable
608               varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
609               DO k=1,klev_src
610                  DO j=1,nbp_lat
611                     DO i=1,nbp_lon
612                        varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k)
613                     END DO
614                  END DO
615               END DO
616             
617               ! Invert latitudes for surface pressure
618               vartmp(:,:) = psurf_glo2D(:,:,imth)
619               DO j=1,nbp_lat
620                  DO i=1,nbp_lon
621                     psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
622                  END DO
623               END DO
624             
625               ! Invert latitudes for the load
626               vartmp(:,:) = load_glo2D(:,:,imth)
627               DO j=1,nbp_lat
628                  DO i=1,nbp_lon
629                     load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
630                  END DO
631               END DO
632            END IF ! invert_lat
633             
634            ! Do zonal mead at poles and distribut at whole first and last latitude
635            DO k=1, klev_src
636               npole=0.  ! North pole, j=1
637               spole=0.  ! South pole, j=nbp_lat       
638               DO i=1,nbp_lon
639                  npole = npole + varyear(i,1,k,imth)
640                  spole = spole + varyear(i,nbp_lat,k,imth)
641               END DO
642               npole = npole/REAL(nbp_lon)
643               spole = spole/REAL(nbp_lon)
644               varyear(:,1,    k,imth) = npole
645               varyear(:,nbp_lat,k,imth) = spole
646            END DO
647         END DO ! imth
648       
649         ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
650         IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1)
651       
652         ! Transform from 2D to 1D field
653         CALL grid2Dto1D_glo(varyear,varyear_glo1D)
654         CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
655         CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
656     
657      ENDIF
658     
659    ELSE
660        ALLOCATE(varyear_glo1D(0,0,0))       
661    END IF ! is_mpi_root .AND. is_omp_root
662
663!$OMP BARRIER
664 
665! 6) Distribute to all processes
666!    Scatter global field(klon_glo) to local process domain(klon)
667!    and distribute klev_src to all processes
668!****************************************************************************************
669
670    ! Distribute klev_src
671    CALL bcast(klev_src)
672
673    ! Allocate and distribute pt_ap and pt_b
674    IF (.NOT. ASSOCIATED(pt_ap)) THEN  ! if pt_ap is allocated also pt_b is allocated
675       ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat=ierr)
676       IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 4',1)
677    END IF
678    CALL bcast(pt_ap)
679    CALL bcast(pt_b)
680
681    ! Allocate space for output pointer variable at local process
682    IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
683    ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
684    ALLOCATE(pt_year_mpi(klon_mpi, klev_src, 12), stat=ierr)
685    IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5',1)
686
687    IF (grid_type==unstructured) THEN
688#ifdef CPP_XIOS
689      IF (is_omp_master) THEN
690        CALL xios_send_field(TRIM(varname)//"_in",varyear)
691        CALL xios_recv_field(TRIM(varname)//"_out",pt_year_mpi)
692        CALL xios_send_field("load_"//TRIM(varname)//"_in",load_glo2D)
693        CALL xios_recv_field("load_"//TRIM(varname)//"_out",load_out_mpi)
694        IF (first) THEN
695          ALLOCATE(psurf_interp(klon_mpi,12))
696          CALL xios_send_field("psurf_aerosol_in",psurf_glo2D)
697          CALL xios_recv_field("psurf_aerosol_out",psurf_interp)
698        ENDIF
699      ENDIF
700      CALL scatter_omp(pt_year_mpi,pt_year)
701      CALL scatter_omp(load_out_mpi,load_out)
702      CALL scatter_omp(psurf_interp,psurf_out)
703      first=.FALSE.
704#endif
705    ELSE
706      ! Scatter global field to local domain at local process
707      CALL scatter(varyear_glo1D, pt_year)
708      CALL scatter(psurf_glo1D, psurf_out)
709      CALL scatter(load_glo1D,  load_out)
710    ENDIF
711! 7) Test for negative values
712!****************************************************************************************
713    IF (MINVAL(pt_year) < 0.) THEN
714       WRITE(lunout,*) 'Warning! Negative values read from file :', fname
715    END IF
716
717  END SUBROUTINE get_aero_fromfile
718
719
720  SUBROUTINE check_err(status,text)
721    USE netcdf
722    USE print_control_mod, ONLY: lunout
723    IMPLICIT NONE
724
725    INTEGER, INTENT (IN) :: status
726    CHARACTER(len=*), INTENT (IN), OPTIONAL :: text
727
728    IF (status /= NF90_NOERR) THEN
729       WRITE(lunout,*) 'Error in get_aero_fromfile, netcdf error code = ',status
730       IF (PRESENT(text)) THEN
731          WRITE(lunout,*) 'Error in get_aero_fromfile : ',text
732       END IF
733       CALL abort_physic('get_aero_fromfile',trim(nf90_strerror(status)),1)
734    END IF
735
736  END SUBROUTINE check_err
737
738
739END MODULE readaerosol_mod
Note: See TracBrowser for help on using the repository browser.