source: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/readaerosol.F90 @ 3655

Last change on this file since 3655 was 1299, checked in by Laurent Fairhead, 15 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

  • tous les FLOAT() sont remplacés par des REAL()
  • tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
  • le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
  • des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper


General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

  • all FLOAT() instructions have been replaced by REAL() instructions
  • all STOP instructions in phylmd have been replaced by calls to abort_gcm
  • the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
  • $Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

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