source: LMDZ5/branches/IPSLCM5A2.1/libf/phylmd/readaerosol.F90 @ 2952

Last change on this file since 2952 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

EM

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