source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/readaerosol.F90 @ 3831

Last change on this file since 3831 was 3831, checked in by ymipsl, 9 years ago

module reorganisation for a cleaner dyn-phys interface
YM

File size: 24.6 KB
Line 
1! $Id: readaerosol.F90 1907 2013-11-26 13:10:46Z lguez $
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, jjp+1) 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
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                                  :: p0           ! Reference pressure value
207    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year      ! Pointer-variabale from file, 12 month, grid : klon,klev_src
208    REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out    ! Surface pression for 12 months
209    REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out     ! Aerosol mass load in each column
210    INTEGER                               :: nbr_tsteps   ! number of month in file read
211
212! Local variables
213    CHARACTER(len=30)     :: fname
214    CHARACTER(len=30)     :: cvar
215    INTEGER               :: ncid, dimid, varid
216    INTEGER               :: imth, i, j, k, ierr
217    REAL                  :: npole, spole
218    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varmth
219    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear       ! Global variable read from file, 12 month
220    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varyear_glo1D !(klon_glo, klev_src, 12)
221    REAL, ALLOCATABLE, DIMENSION(:)       :: varktmp
222
223    REAL, DIMENSION(nbp_lon,nbp_lat,12)   :: psurf_glo2D   ! Surface pression for 12 months on dynamics global grid
224    REAL, DIMENSION(klon_glo,12)          :: psurf_glo1D   ! -"- on physical global grid
225    REAL, DIMENSION(nbp_lon,nbp_lat,12)   :: load_glo2D    ! Load for 12 months on dynamics global grid
226    REAL, DIMENSION(klon_glo,12)          :: load_glo1D    ! -"- on physical global grid
227    REAL, DIMENSION(nbp_lon,nbp_lat)      :: vartmp
228    REAL, DIMENSION(nbp_lon)              :: lon_src              ! longitudes in file
229    REAL, DIMENSION(nbp_lat)              :: lat_src, lat_src_inv ! latitudes in file
230    LOGICAL                               :: new_file             ! true if new file format detected
231    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
232
233
234    ! Deallocate pointers
235    IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
236    IF (ASSOCIATED(pt_b))  DEALLOCATE(pt_b)
237
238    IF (is_mpi_root .AND. is_omp_root) THEN
239
240! 1) Open file
241!****************************************************************************************
242! Add suffix to filename
243       fname = trim(filename)//cyr//'.nc'
244 
245       WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname)
246       CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(varname) )
247
248! Test for equal longitudes and latitudes in file and model
249!****************************************************************************************
250       ! Read and test longitudes
251       CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
252       CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
253       
254       IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
255          WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
256          WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
257          WRITE(lunout,*) 'longitudes in model :', io_lon
258         
259          CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
260       END IF
261
262       ! Read and test latitudes
263       CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
264       CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
265
266       ! Invert source latitudes
267       DO j = 1, nbp_lat
268          lat_src_inv(j) = lat_src(nbp_lat +1 -j)
269       END DO
270
271       IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
272          ! Latitudes are the same
273          invert_lat=.FALSE.
274       ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
275          ! Inverted source latitudes correspond to model latitudes
276          WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
277          invert_lat=.TRUE.
278       ELSE
279          WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
280          WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src     
281          WRITE(lunout,*) 'latitudes in model :', io_lat
282          CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
283       END IF
284
285
286! 2) Check if old or new file is avalabale.
287!    New type of file should contain the dimension 'lev'
288!    Old type of file should contain the dimension 'PRESNIVS'
289!****************************************************************************************
290       ierr = nf90_inq_dimid(ncid, 'lev', dimid)
291       IF (ierr /= NF90_NOERR) THEN
292          ! Coordinate axe lev not found. Check for presnivs.
293          ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid)
294          IF (ierr /= NF90_NOERR) THEN
295             ! Dimension PRESNIVS not found either
296             CALL abort_physic('get_aero_fromfile', 'dimension lev or presnivs not in file',1)
297          ELSE
298             ! Old file found
299             new_file=.FALSE.
300             WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done'
301          END IF
302       ELSE
303          ! New file found
304          new_file=.TRUE.
305          WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
306       END IF
307       
308! 2) Find vertical dimension klev_src
309!****************************************************************************************
310       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
311       
312     ! Allocate variables depending on the number of vertical levels
313       ALLOCATE(varmth(nbp_lon, nbp_lat, klev_src), varyear(nbp_lon, nbp_lat, klev_src, 12), stat=ierr)
314       IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 1',1)
315
316       ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr)
317       IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 2',1)
318
319! 3) Read all variables from file
320!    There is 2 options for the file structure :
321!    new_file=TRUE  : read varyear, ps, pt_ap and pt_b
322!    new_file=FALSE : read varyear month by month
323!****************************************************************************************
324
325       IF (new_file) THEN
326! ++) Check number of month in file opened
327!**************************************************************************************************
328       ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
329       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps),"pb inq dim TIME" )
330!       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
331       IF (nbr_tsteps /= 12 ) THEN
332         CALL abort_physic('get_aero_fromfile', &
333         'not the right number of months in aerosol file read (should be 12 for the moment)',1)
334       ENDIF
335
336! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
337!****************************************************************************************
338          ! Get variable id
339          CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) )
340         
341          ! Get the variable
342          CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)),"pb get var "//TRIM(varname) )
343         
344! ++) Read surface pression, 12 month in one variable
345!****************************************************************************************
346          ! Get variable id
347          CALL check_err( nf90_inq_varid(ncid, "ps", varid),"pb inq var ps" )
348          ! Get the variable
349          CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D),"pb get var ps" )
350         
351! ++) Read mass load, 12 month in one variable
352!****************************************************************************************
353          ! Get variable id
354          CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname))
355          ! Get the variable
356          CALL check_err( nf90_get_var(ncid, varid, load_glo2D),"pb get var load_"//TRIM(varname) )
357         
358! ++) Read ap
359!****************************************************************************************
360          ! Get variable id
361          CALL check_err( nf90_inq_varid(ncid, "ap", varid),"pb inq var ap" )
362          ! Get the variable
363          CALL check_err( nf90_get_var(ncid, varid, pt_ap),"pb get var ap" )
364
365! ++) Read b
366!****************************************************************************************
367          ! Get variable id
368          CALL check_err( nf90_inq_varid(ncid, "b", varid),"pb inq var b" )
369          ! Get the variable
370          CALL check_err( nf90_get_var(ncid, varid, pt_b),"pb get var b" )
371
372! ++) Read p0 : reference pressure
373!****************************************************************************************
374          ! Get variable id
375          CALL check_err( nf90_inq_varid(ncid, "p0", varid),"pb inq var p0" )
376          ! Get the variable
377          CALL check_err( nf90_get_var(ncid, varid, p0),"pb get var p0" )
378         
379
380       ELSE  ! old file
381
382! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
383!****************************************************************************************
384          DO imth=1, 12
385             IF (imth.EQ.1) THEN
386                cvar=TRIM(varname)//'JAN'
387             ELSE IF (imth.EQ.2) THEN
388                cvar=TRIM(varname)//'FEB'
389             ELSE IF (imth.EQ.3) THEN
390                cvar=TRIM(varname)//'MAR'
391             ELSE IF (imth.EQ.4) THEN
392                cvar=TRIM(varname)//'APR'
393             ELSE IF (imth.EQ.5) THEN
394                cvar=TRIM(varname)//'MAY'
395             ELSE IF (imth.EQ.6) THEN
396                cvar=TRIM(varname)//'JUN'
397             ELSE IF (imth.EQ.7) THEN
398                cvar=TRIM(varname)//'JUL'
399             ELSE IF (imth.EQ.8) THEN
400                cvar=TRIM(varname)//'AUG'
401             ELSE IF (imth.EQ.9) THEN
402                cvar=TRIM(varname)//'SEP'
403             ELSE IF (imth.EQ.10) THEN
404                cvar=TRIM(varname)//'OCT'
405             ELSE IF (imth.EQ.11) THEN
406                cvar=TRIM(varname)//'NOV'
407             ELSE IF (imth.EQ.12) THEN
408                cvar=TRIM(varname)//'DEC'
409             END IF
410             
411             ! Get variable id
412             CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid),"pb inq var "//TRIM(cvar) )
413             
414             ! Get the variable
415             CALL check_err( nf90_get_var(ncid, varid, varmth),"pb get var "//TRIM(cvar) )
416             
417             ! Store in variable for the whole year
418             varyear(:,:,:,imth)=varmth(:,:,:)
419             
420          END DO
421         
422          ! Putting dummy
423          psurf_glo2D(:,:,:) = not_valid
424          load_glo2D(:,:,:)  = not_valid
425          pt_ap(:) = not_valid
426          pt_b(:)  = not_valid
427
428       END IF
429
430! 4) Close file 
431!****************************************************************************************
432       CALL check_err( nf90_close(ncid),"pb in close" )
433     
434
435! 5) Transform the global field from 2D(nbp_lon, jjp+1) to 1D(klon_glo)
436!****************************************************************************************
437! Test if vertical levels have to be inversed
438
439       IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN
440!          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inverted'
441!          WRITE(lunout,*) 'before pt_ap = ', pt_ap
442!          WRITE(lunout,*) 'before pt_b = ', pt_b
443         
444          ! Inverse vertical levels for varyear
445          DO imth=1, 12
446             varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
447             DO k=1, klev_src
448                DO j=1, nbp_lat
449                   DO i=1,nbp_lon
450                      varyear(i,j,k,imth) = varmth(i,j,klev_src+1-k)
451                   END DO
452                END DO
453             END DO
454          END DO
455           
456          ! Inverte vertical axes for pt_ap and pt_b
457          varktmp(:) = pt_ap(:)
458          DO k=1, klev_src
459             pt_ap(k) = varktmp(klev_src+1-k)
460          END DO
461
462          varktmp(:) = pt_b(:)
463          DO k=1, klev_src
464             pt_b(k) = varktmp(klev_src+1-k)
465          END DO
466          WRITE(lunout,*) 'after pt_ap = ', pt_ap
467          WRITE(lunout,*) 'after pt_b = ', pt_b
468
469       ELSE
470          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
471          WRITE(lunout,*) 'pt_ap = ', pt_ap
472          WRITE(lunout,*) 'pt_b = ', pt_b
473       END IF
474
475!     - Invert latitudes if necessary
476       DO imth=1, 12
477          IF (invert_lat) THEN
478
479             ! Invert latitudes for the variable
480             varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
481             DO k=1,klev_src
482                DO j=1,nbp_lat
483                   DO i=1,nbp_lon
484                      varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k)
485                   END DO
486                END DO
487             END DO
488             
489             ! Invert latitudes for surface pressure
490             vartmp(:,:) = psurf_glo2D(:,:,imth)
491             DO j=1, nbp_lat
492                DO i=1,nbp_lon
493                   psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
494                END DO
495             END DO
496             
497             ! Invert latitudes for the load
498             vartmp(:,:) = load_glo2D(:,:,imth)
499             DO j=1, nbp_lat
500                DO i=1,nbp_lon
501                   load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
502                END DO
503             END DO
504          END IF ! invert_lat
505             
506          ! Do zonal mead at poles and distribut at whole first and last latitude
507          DO k=1, klev_src
508             npole=0.  ! North pole, j=1
509             spole=0.  ! South pole, j=nbp_lat         
510             DO i=1,nbp_lon
511                npole = npole + varyear(i,1,k,imth)
512                spole = spole + varyear(i,nbp_lat,k,imth)
513             END DO
514             npole = npole/REAL(nbp_lon)
515             spole = spole/REAL(nbp_lon)
516             varyear(:,1,    k,imth) = npole
517             varyear(:,nbp_lat,k,imth) = spole
518          END DO
519       END DO ! imth
520       
521       ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
522       IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1)
523       
524       ! Transform from 2D to 1D field
525       CALL grid2Dto1D_glo(varyear,varyear_glo1D)
526       CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
527       CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
528
529    ELSE
530      ALLOCATE(varyear_glo1D(0,0,0))       
531    END IF ! is_mpi_root .AND. is_omp_root
532
533!$OMP BARRIER
534 
535! 6) Distribute to all processes
536!    Scatter global field(klon_glo) to local process domain(klon)
537!    and distribute klev_src to all processes
538!****************************************************************************************
539
540    ! Distribute klev_src
541    CALL bcast(klev_src)
542
543    ! Allocate and distribute pt_ap and pt_b
544    IF (.NOT. ASSOCIATED(pt_ap)) THEN  ! if pt_ap is allocated also pt_b is allocated
545       ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat=ierr)
546       IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 4',1)
547    END IF
548    CALL bcast(pt_ap)
549    CALL bcast(pt_b)
550
551    ! Allocate space for output pointer variable at local process
552    IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
553    ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
554    IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5',1)
555
556    ! Scatter global field to local domain at local process
557    CALL scatter(varyear_glo1D, pt_year)
558    CALL scatter(psurf_glo1D, psurf_out)
559    CALL scatter(load_glo1D,  load_out)
560
561! 7) Test for negative values
562!****************************************************************************************
563    IF (MINVAL(pt_year) < 0.) THEN
564       WRITE(lunout,*) 'Warning! Negative values read from file :', fname
565    END IF
566
567  END SUBROUTINE get_aero_fromfile
568
569
570  SUBROUTINE check_err(status,text)
571    USE netcdf
572    USE print_control_mod, ONLY: lunout
573    IMPLICIT NONE
574
575    INTEGER, INTENT (IN) :: status
576    CHARACTER(len=*), INTENT (IN), OPTIONAL :: text
577
578    IF (status /= NF90_NOERR) THEN
579       WRITE(lunout,*) 'Error in get_aero_fromfile, netcdf error code = ',status
580       IF (PRESENT(text)) THEN
581          WRITE(lunout,*) 'Error in get_aero_fromfile : ',text
582       END IF
583       CALL abort_physic('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.