source: LMDZ5/branches/IPSLCM6.0.10/libf/phylmd/readaerosol.F90 @ 3312

Last change on this file since 3312 was 2845, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2838:2842 into testing branch

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