source: LMDZ5/branches/LMDZ5V1.0-dev/libf/phylmd/readaerosol.F90 @ 3619

Last change on this file since 3619 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

  • IPSLCM5A configuration (comparison between trunk and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent
  • MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

  • IPSLCM5A configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent
  • MH07 configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files


Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

  • pour une config IPSLCM5A (comparaison tronc/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux (à part sur RN et Pb)
    • fichiers histoire égaux
  • pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux
    • fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

  • pour une config IPSLCM5A:
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc
  • pour une config MH07
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.6 KB
Line 
1! $Id: readaerosol.F90 1403 2010-07-01 09:02:53Z fhourdin $
2!
3MODULE readaerosol_mod
4
5  REAL, SAVE :: not_valid=-333.
6
7CONTAINS
8
9SUBROUTINE readaerosol(name_aero, type, 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  ! correspond to aer_type in clesphys.h
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, 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, 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, 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, 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, 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, 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, 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_gcm('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
163     CALL abort_gcm('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, 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(iim, 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
194    IMPLICIT NONE
195     
196    INCLUDE "dimensions.h"     
197    INCLUDE "iniprint.h"
198
199! Input argumets
200    CHARACTER(len=7), INTENT(IN)          :: varname
201    CHARACTER(len=4), INTENT(IN)          :: cyr
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=8)      :: filename='aerosols'
216    CHARACTER(len=30)     :: cvar
217    INTEGER               :: ncid, dimid, varid
218    INTEGER               :: imth, i, j, k, ierr
219    REAL                  :: npole, spole
220    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varmth
221    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear       ! Global variable read from file, 12 month
222    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varyear_glo1D !(klon_glo, klev_src, 12)
223    REAL, ALLOCATABLE, DIMENSION(:)       :: varktmp
224
225    REAL, DIMENSION(iim,jjm+1,12)         :: psurf_glo2D   ! Surface pression for 12 months on dynamics global grid
226    REAL, DIMENSION(klon_glo,12)          :: psurf_glo1D   ! -"- on physical global grid
227    REAL, DIMENSION(iim,jjm+1,12)         :: load_glo2D    ! Load for 12 months on dynamics global grid
228    REAL, DIMENSION(klon_glo,12)          :: load_glo1D    ! -"- on physical global grid
229    REAL, DIMENSION(iim,jjm+1)            :: vartmp
230    REAL, DIMENSION(iim)                  :: lon_src              ! longitudes in file
231    REAL, DIMENSION(jjm+1)                :: lat_src, lat_src_inv ! latitudes in file
232    LOGICAL                               :: new_file             ! true if new file format detected
233    LOGICAL                               :: invert_lat           ! true if the field has to be inverted for latitudes
234
235
236    ! Deallocate pointers
237    IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
238    IF (ASSOCIATED(pt_b))  DEALLOCATE(pt_b)
239
240    IF (is_mpi_root .AND. is_omp_root) THEN
241
242! 1) Open file
243!****************************************************************************************
244       fname = filename//cyr//'.nc'
245 
246       WRITE(lunout,*) 'reading ', TRIM(fname)
247       CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid) )
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) )
253       CALL check_err( nf90_get_var(ncid, varid, lon_src(:)) )
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_gcm('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) )
265       CALL check_err( nf90_get_var(ncid, varid, lat_src(:)) )
266
267       ! Invert source latitudes
268       DO j = 1, jjm+1
269          lat_src_inv(j) = lat_src(jjm+1 +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_gcm('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
284       END IF
285
286! 1.5) Check number of month in file opened
287!
288!**************************************************************************************************
289       ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
290       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps) )
291!       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
292       IF (nbr_tsteps /= 12 ) THEN
293         CALL abort_gcm('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)',1)
294       ENDIF
295
296
297! 2) Check if old or new file is avalabale.
298!    New type of file should contain the dimension 'lev'
299!    Old type of file should contain the dimension 'PRESNIVS'
300!****************************************************************************************
301       ierr = nf90_inq_dimid(ncid, 'lev', dimid)
302       IF (ierr /= NF90_NOERR) THEN
303          ! Coordinate axe lev not found. Check for presnivs.
304          ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid)
305          IF (ierr /= NF90_NOERR) THEN
306             ! Dimension PRESNIVS not found either
307             CALL abort_gcm('get_aero_fromfile', 'dimension lev or presnivs not in file',1)
308          ELSE
309             ! Old file found
310             new_file=.FALSE.
311             WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done'
312          END IF
313       ELSE
314          ! New file found
315          new_file=.TRUE.
316          WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
317       END IF
318       
319! 2) Find vertical dimension klev_src
320!****************************************************************************************
321       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src) )
322       
323     ! Allocate variables depending on the number of vertical levels
324       ALLOCATE(varmth(iim, jjm+1, klev_src), varyear(iim, jjm+1, klev_src, 12), stat=ierr)
325       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 1',1)
326
327       ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr)
328       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 2',1)
329
330! 3) Read all variables from file
331!    There is 2 options for the file structure :
332!    new_file=TRUE  : read varyear, ps, pt_ap and pt_b
333!    new_file=FALSE : read varyear month by month
334!****************************************************************************************
335
336       IF (new_file) THEN
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) )
342         
343          ! Get the variable
344          CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)) )
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) )
350          ! Get the variable
351          CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D) )
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) )
357          ! Get the variable
358          CALL check_err( nf90_get_var(ncid, varid, load_glo2D) )
359         
360! ++) Read ap
361!****************************************************************************************
362          ! Get variable id
363          CALL check_err( nf90_inq_varid(ncid, "ap", varid) )
364          ! Get the variable
365          CALL check_err( nf90_get_var(ncid, varid, pt_ap) )
366
367! ++) Read b
368!****************************************************************************************
369          ! Get variable id
370          CALL check_err( nf90_inq_varid(ncid, "b", varid) )
371          ! Get the variable
372          CALL check_err( nf90_get_var(ncid, varid, pt_b) )
373
374! ++) Read p0 : reference pressure
375!****************************************************************************************
376          ! Get variable id
377          CALL check_err( nf90_inq_varid(ncid, "p0", varid) )
378          ! Get the variable
379          CALL check_err( nf90_get_var(ncid, varid, 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) )
415             
416             ! Get the variable
417             CALL check_err( nf90_get_var(ncid, varid, varmth) )
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) )
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)
573    USE netcdf
574    IMPLICIT NONE
575
576    INCLUDE "iniprint.h"
577    INTEGER, INTENT (IN) :: status
578
579    IF (status /= NF90_NOERR) THEN
580       WRITE(lunout,*) 'Error in get_aero_fromfile ',status
581       CALL abort_gcm('get_aero_fromfile',trim(nf90_strerror(status)),1)
582    END IF
583
584  END SUBROUTINE check_err
585
586
587END MODULE readaerosol_mod
Note: See TracBrowser for help on using the repository browser.