source: LMDZ6/branches/blowing_snow/libf/phylmd/readaerosol_mod.F90

Last change on this file was 4103, checked in by Laurent Fairhead, 3 years ago

Inclusion of some corrections and optimisations for XIOS done by
Arnaud Durocher during his TGCC mission.
Included here are r3703, r3704, r3750, r3751, r3752 from his
LMDZ6/branches/Optimisation_LMDZ branch

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