source: LMDZ4/branches/LMDZ4_AR5/libf/phylmd/readaerosol.F90 @ 1483

Last change on this file since 1483 was 1483, checked in by jghattas, 13 years ago
  • Added option aer_type=mix1 and mix2

mix1: as aer_type=scenario for SO4 and aer_type=annuel for the other aerosols.
mix2: as aer_type=scenario for SO4 and aer_type=preind for the other aerosols.
For these 2 options 2 files are read : so4.runxxxx.nc for SO4 and aerosolsxxxx.nc the other aerosols.

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