source: lmdz_wrf/WRFV3/lmdz/readaerosol.F90 @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 24.5 KB
Line 
1! $Id: readaerosol.F90 1665 2012-10-09 13:35:26Z 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
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) - REAL(iyr_in-iyr1)/REAL(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) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * &
148                      (psurf(i,it) - psurf2(i,it))
149
150                 load(i,it) = &
151                      load(i,it) - REAL(iyr_in-iyr1)/REAL(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), "pb open "//trim(varname) )
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),"pb inq lon" )
255       CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
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),"pb inq lat" )
267       CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
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),"pb inq dim for PRESNIVS or lev" )
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! ++) Check number of month in file opened
330!**************************************************************************************************
331       ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
332       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps),"pb inq dim TIME" )
333!       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
334       IF (nbr_tsteps /= 12 ) THEN
335         CALL abort_gcm('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)',1)
336       ENDIF
337
338! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
339!****************************************************************************************
340          ! Get variable id
341          CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) )
342         
343          ! Get the variable
344          CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)),"pb get var "//TRIM(varname) )
345         
346! ++) Read surface pression, 12 month in one variable
347!****************************************************************************************
348          ! Get variable id
349          CALL check_err( nf90_inq_varid(ncid, "ps", varid),"pb inq var ps" )
350          ! Get the variable
351          CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D),"pb get var ps" )
352         
353! ++) Read mass load, 12 month in one variable
354!****************************************************************************************
355          ! Get variable id
356          CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname))
357          ! Get the variable
358          CALL check_err( nf90_get_var(ncid, varid, load_glo2D),"pb get var load_"//TRIM(varname) )
359         
360! ++) Read ap
361!****************************************************************************************
362          ! Get variable id
363          CALL check_err( nf90_inq_varid(ncid, "ap", varid),"pb inq var ap" )
364          ! Get the variable
365          CALL check_err( nf90_get_var(ncid, varid, pt_ap),"pb get var ap" )
366
367! ++) Read b
368!****************************************************************************************
369          ! Get variable id
370          CALL check_err( nf90_inq_varid(ncid, "b", varid),"pb inq var b" )
371          ! Get the variable
372          CALL check_err( nf90_get_var(ncid, varid, pt_b),"pb get var b" )
373
374! ++) Read p0 : reference pressure
375!****************************************************************************************
376          ! Get variable id
377          CALL check_err( nf90_inq_varid(ncid, "p0", varid),"pb inq var p0" )
378          ! Get the variable
379          CALL check_err( nf90_get_var(ncid, varid, p0),"pb get var p0" )
380         
381
382       ELSE  ! old file
383
384! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
385!****************************************************************************************
386          DO imth=1, 12
387             IF (imth.EQ.1) THEN
388                cvar=TRIM(varname)//'JAN'
389             ELSE IF (imth.EQ.2) THEN
390                cvar=TRIM(varname)//'FEB'
391             ELSE IF (imth.EQ.3) THEN
392                cvar=TRIM(varname)//'MAR'
393             ELSE IF (imth.EQ.4) THEN
394                cvar=TRIM(varname)//'APR'
395             ELSE IF (imth.EQ.5) THEN
396                cvar=TRIM(varname)//'MAY'
397             ELSE IF (imth.EQ.6) THEN
398                cvar=TRIM(varname)//'JUN'
399             ELSE IF (imth.EQ.7) THEN
400                cvar=TRIM(varname)//'JUL'
401             ELSE IF (imth.EQ.8) THEN
402                cvar=TRIM(varname)//'AUG'
403             ELSE IF (imth.EQ.9) THEN
404                cvar=TRIM(varname)//'SEP'
405             ELSE IF (imth.EQ.10) THEN
406                cvar=TRIM(varname)//'OCT'
407             ELSE IF (imth.EQ.11) THEN
408                cvar=TRIM(varname)//'NOV'
409             ELSE IF (imth.EQ.12) THEN
410                cvar=TRIM(varname)//'DEC'
411             END IF
412             
413             ! Get variable id
414             CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid),"pb inq var "//TRIM(cvar) )
415             
416             ! Get the variable
417             CALL check_err( nf90_get_var(ncid, varid, varmth),"pb get var "//TRIM(cvar) )
418             
419             ! Store in variable for the whole year
420             varyear(:,:,:,imth)=varmth(:,:,:)
421             
422          END DO
423         
424          ! Putting dummy
425          psurf_glo2D(:,:,:) = not_valid
426          load_glo2D(:,:,:)  = not_valid
427          pt_ap(:) = not_valid
428          pt_b(:)  = not_valid
429
430       END IF
431
432! 4) Close file 
433!****************************************************************************************
434       CALL check_err( nf90_close(ncid),"pb in close" )
435     
436
437! 5) Transform the global field from 2D(iim, jjp+1) to 1D(klon_glo)
438!****************************************************************************************
439! Test if vertical levels have to be inversed
440
441       IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN
442!          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inverted'
443!          WRITE(lunout,*) 'before pt_ap = ', pt_ap
444!          WRITE(lunout,*) 'before pt_b = ', pt_b
445         
446          ! Inverse vertical levels for varyear
447          DO imth=1, 12
448             varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
449             DO k=1, klev_src
450                DO j=1, jjm+1
451                   DO i=1,iim
452                      varyear(i,j,k,imth) = varmth(i,j,klev_src+1-k)
453                   END DO
454                END DO
455             END DO
456          END DO
457           
458          ! Inverte vertical axes for pt_ap and pt_b
459          varktmp(:) = pt_ap(:)
460          DO k=1, klev_src
461             pt_ap(k) = varktmp(klev_src+1-k)
462          END DO
463
464          varktmp(:) = pt_b(:)
465          DO k=1, klev_src
466             pt_b(k) = varktmp(klev_src+1-k)
467          END DO
468          WRITE(lunout,*) 'after pt_ap = ', pt_ap
469          WRITE(lunout,*) 'after pt_b = ', pt_b
470
471       ELSE
472          WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
473          WRITE(lunout,*) 'pt_ap = ', pt_ap
474          WRITE(lunout,*) 'pt_b = ', pt_b
475       END IF
476
477!     - Invert latitudes if necessary
478       DO imth=1, 12
479          IF (invert_lat) THEN
480
481             ! Invert latitudes for the variable
482             varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
483             DO k=1,klev_src
484                DO j=1,jjm+1
485                   DO i=1,iim
486                      varyear(i,j,k,imth) = varmth(i,jjm+1+1-j,k)
487                   END DO
488                END DO
489             END DO
490             
491             ! Invert latitudes for surface pressure
492             vartmp(:,:) = psurf_glo2D(:,:,imth)
493             DO j=1, jjm+1
494                DO i=1,iim
495                   psurf_glo2D(i,j,imth)= vartmp(i,jjm+1+1-j)
496                END DO
497             END DO
498             
499             ! Invert latitudes for the load
500             vartmp(:,:) = load_glo2D(:,:,imth)
501             DO j=1, jjm+1
502                DO i=1,iim
503                   load_glo2D(i,j,imth)= vartmp(i,jjm+1+1-j)
504                END DO
505             END DO
506          END IF ! invert_lat
507             
508          ! Do zonal mead at poles and distribut at whole first and last latitude
509          DO k=1, klev_src
510             npole=0.  ! North pole, j=1
511             spole=0.  ! South pole, j=jjm+1         
512             DO i=1,iim
513                npole = npole + varyear(i,1,k,imth)
514                spole = spole + varyear(i,jjm+1,k,imth)
515             END DO
516             npole = npole/REAL(iim)
517             spole = spole/REAL(iim)
518             varyear(:,1,    k,imth) = npole
519             varyear(:,jjm+1,k,imth) = spole
520          END DO
521       END DO ! imth
522       
523       ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
524       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 3',1)
525       
526       ! Transform from 2D to 1D field
527       CALL grid2Dto1D_glo(varyear,varyear_glo1D)
528       CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
529       CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
530
531    ELSE
532      ALLOCATE(varyear_glo1D(0,0,0))       
533    END IF ! is_mpi_root .AND. is_omp_root
534
535!$OMP BARRIER
536 
537! 6) Distribute to all processes
538!    Scatter global field(klon_glo) to local process domain(klon)
539!    and distribute klev_src to all processes
540!****************************************************************************************
541
542    ! Distribute klev_src
543    CALL bcast(klev_src)
544
545    ! Allocate and distribute pt_ap and pt_b
546    IF (.NOT. ASSOCIATED(pt_ap)) THEN  ! if pt_ap is allocated also pt_b is allocated
547       ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat=ierr)
548       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 4',1)
549    END IF
550    CALL bcast(pt_ap)
551    CALL bcast(pt_b)
552
553    ! Allocate space for output pointer variable at local process
554    IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
555    ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
556    IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5',1)
557
558    ! Scatter global field to local domain at local process
559    CALL scatter(varyear_glo1D, pt_year)
560    CALL scatter(psurf_glo1D, psurf_out)
561    CALL scatter(load_glo1D,  load_out)
562
563! 7) Test for negative values
564!****************************************************************************************
565    IF (MINVAL(pt_year) < 0.) THEN
566       WRITE(lunout,*) 'Warning! Negative values read from file :', fname
567    END IF
568
569  END SUBROUTINE get_aero_fromfile
570
571
572  SUBROUTINE check_err(status,text)
573    USE netcdf
574    IMPLICIT NONE
575
576    INCLUDE "iniprint.h"
577    INTEGER, INTENT (IN) :: status
578    CHARACTER(len=*), INTENT (IN), OPTIONAL :: text
579
580    IF (status /= NF90_NOERR) THEN
581       WRITE(lunout,*) 'Error in get_aero_fromfile, netcdf error code = ',status
582       IF (PRESENT(text)) THEN
583          WRITE(lunout,*) 'Error in get_aero_fromfile : ',text
584       END IF
585       CALL abort_gcm('get_aero_fromfile',trim(nf90_strerror(status)),1)
586    END IF
587
588  END SUBROUTINE check_err
589
590
591END MODULE readaerosol_mod
Note: See TracBrowser for help on using the repository browser.