Index: LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90
===================================================================
--- LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90	(revision 1249)
+++ LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90	(revision 1265)
@@ -7,5 +7,5 @@
 CONTAINS
 
-SUBROUTINE readaerosol(name_aero, type, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
+SUBROUTINE readaerosol(name_aero, type, iyr_in, klev_src, pt_ap, pt_b, pt_out, pt_psurf, pt_load, nbr_tsteps)
 
 !****************************************************************************************
@@ -34,16 +34,17 @@
   REAL, POINTER, DIMENSION(:)     :: pt_ap        ! Pointer for describing the vertical levels      
   REAL, POINTER, DIMENSION(:)     :: pt_b         ! Pointer for describing the vertical levels      
-  REAL, POINTER, DIMENSION(:,:,:) :: pt_out       ! The massvar distributions, DIMENSION(klon, klev_src, 12)
-  REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf  ! Surface pression for 12 months
-  REAL, DIMENSION(klon,12), INTENT(OUT) :: load   ! Aerosol mass load in each column for 12 months
+  REAL, POINTER, DIMENSION(:,:,:) :: pt_out       ! The massvar distributions, DIMENSION(klon,klev_src, nbr_tsteps)
+  REAL, POINTER, DIMENSION(:,:)   :: pt_psurf     ! The massvar distributions, DIMENSION(klon, nbr_tsteps)
+  REAL, POINTER, DIMENSION(:,:)   :: pt_load      ! The massvar distributions, DIMENSION(klon, nbr_tsteps)
+  INTEGER                         :: nbr_tsteps   ! number of timesteps in read file (12 or 14)
 
   ! Local variables
-  CHARACTER(len=4)                :: cyear
-  REAL, POINTER, DIMENSION(:,:,:) :: pt_2
-  REAL, DIMENSION(klon,12)        :: psurf2, load2
-  REAL                            :: p0           ! Reference pressure
-  INTEGER                         :: iyr1, iyr2, klev_src2
-  INTEGER                         :: it, k, i
-  LOGICAL, PARAMETER              :: lonlyone=.FALSE.
+  CHARACTER(len=4)                  :: cyear
+  REAL, POINTER, DIMENSION(:,:,:)   :: pt_2
+  REAL, POINTER, DIMENSION(:,:)     :: psurf2, load2
+  REAL                              :: p0           ! Reference pressure
+  INTEGER                           :: iyr1, iyr2, klev_src2
+  INTEGER                           :: it, k, i
+  LOGICAL, PARAMETER                :: lonlyone=.FALSE.
 
 !****************************************************************************************
@@ -56,7 +57,7 @@
 !****************************************************************************************
      cyear='1980'
-     ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month 
-     ! pt_out has dimensions (klon, klev_src, 12)
-     CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
+     ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 
+     ! pt_out has dimensions (klon, klev_src, nbr_tsteps)
+     CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
      
 
@@ -65,7 +66,17 @@
 !****************************************************************************************     
      cyear='.nat'
-     ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month 
-     ! pt_out has dimensions (klon, klev_src, 12)
-     CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
+     ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 
+     ! pt_out has dimensions (klon, klev_src, nbr_tsteps)
+     CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
+
+  ELSE IF (type == 'annuel') THEN
+! Read and return data from scenario annual files
+!****************************************************************************************     
+     WRITE(cyear,'(I4)') iyr_in
+     WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,'   ',cyear
+     ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 
+     ! pt_out has dimensions (klon, klev_src, nbr_tsteps)
+     CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
+     
      
   ELSE IF (type == 'scenario') THEN
@@ -75,14 +86,14 @@
         cyear='.nat'
         WRITE(lunout,*) 'get_aero 1 iyr_in=', iyr_in,'   ',cyear
-        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month 
-        ! pt_out has dimensions (klon, klev_src, 12)
-        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
+        ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 
+        ! pt_out has dimensions (klon, klev_src, nbr_tsteps)
+        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
         
      ELSE IF (iyr_in .GE. 2100) THEN
         cyear='2100'
         WRITE(lunout,*) 'get_aero 2 iyr_in=', iyr_in,'   ',cyear
-        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month 
-        ! pt_out has dimensions (klon, klev_src, 12)
-        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
+        ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 
+        ! pt_out has dimensions (klon, klev_src, nbr_tsteps)
+        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
         
      ELSE
@@ -102,7 +113,7 @@
         WRITE(cyear,'(I4)') iyr1
         WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,'   ',cyear
-        ! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month 
-        ! pt_out has dimensions (klon, klev_src, 12)
-        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
+        ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month 
+        ! pt_out has dimensions (klon, klev_src, nbr_tsteps)
+        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
         
         ! If to read two decades:
@@ -114,7 +125,7 @@
            
            NULLIFY(pt_2)
-           ! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month 
-           ! pt_2 has dimensions (klon, klev_src, 12)
-           CALL get_aero_fromfile(name_aero, cyear, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2)
+           ! get_aero_fromfile returns pt_2 allocated and initialized with data for nbr_tsteps month 
+           ! pt_2 has dimensions (klon, klev_src, nbr_tsteps)
+           CALL get_aero_fromfile(name_aero, cyear, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2, nbr_tsteps)
            ! Test for same number of vertical levels
            IF (klev_src /= klev_src2) THEN
@@ -124,5 +135,5 @@
            
            ! Linare interpolate to the actual year:
-           DO it=1,12
+           DO it=1,nbr_tsteps
               DO k=1,klev_src
                  DO i = 1, klon
@@ -134,11 +145,11 @@
 
               DO i = 1, klon
-                 psurf(i,it) = &
-                      psurf(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
-                      (psurf(i,it) - psurf2(i,it))
-
-                 load(i,it) = &
-                      load(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
-                      (load(i,it) - load2(i,it))
+                 pt_psurf(i,it) = &
+                      pt_psurf(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
+                      (pt_psurf(i,it) - psurf2(i,it))
+
+                 pt_load(i,it) = &
+                      pt_load(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
+                      (pt_load(i,it) - load2(i,it))
               END DO
            END DO
@@ -146,4 +157,6 @@
            ! Deallocate pt_2 no more needed
            DEALLOCATE(pt_2)
+           DEALLOCATE(psurf2)
+           DEALLOCATE(load2)
            
         END IF ! lonlyone
@@ -159,7 +172,7 @@
 
 
-  SUBROUTINE get_aero_fromfile(varname, cyr, klev_src, pt_ap, pt_b, p0, pt_year, psurf_out, load_out)
-!****************************************************************************************
-! Read 12 month aerosol from file and distribute to local process on physical grid. 
+  SUBROUTINE get_aero_fromfile(varname, cyr, klev_src, pt_ap, pt_b, p0, pt_year, pt_psurf_out, pt_load_out, nbr_tsteps)
+!****************************************************************************************
+! Read nbr_tsteps month aerosol from file and distribute to local process on physical grid. 
 ! Vertical levels, klev_src, may differ from model levels if new file format.
 !
@@ -197,7 +210,8 @@
     REAL, POINTER, DIMENSION(:)           :: pt_b         ! Pointer for describing the vertical levels      
     REAL                                  :: p0           ! Reference pressure value
-    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year      ! Pointer-variabale from file, 12 month, grid : klon,klev_src
-    REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out    ! Surface pression for 12 months
-    REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out     ! Aerosol mass load in each column
+    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year      ! Pointer-variabale from file, nbr_tsteps month, grid : klon,klev_src
+    REAL, POINTER, DIMENSION(:,:)         :: pt_psurf_out    ! Surface pression for nbr_tsteps months
+    REAL, POINTER, DIMENSION(:,:)         :: pt_load_out     ! Aerosol mass load in each column
+    INTEGER                               :: nbr_tsteps ! number of month in file read
 
 ! Local variables
@@ -209,12 +223,12 @@
     REAL                  :: npole, spole
     REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varmth
-    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear       ! Global variable read from file, 12 month
-    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varyear_glo1D !(klon_glo, klev_src, 12)
+    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear       ! Global variable read from file, nbr_tsteps month
+    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varyear_glo1D !(klon_glo, klev_src, nbr_tsteps)
     REAL, ALLOCATABLE, DIMENSION(:)       :: varktmp
 
-    REAL, DIMENSION(iim,jjm+1,12)         :: psurf_glo2D   ! Surface pression for 12 months on dynamics global grid
-    REAL, DIMENSION(klon_glo,12)          :: psurf_glo1D   ! -"- on physical global grid
-    REAL, DIMENSION(iim,jjm+1,12)         :: load_glo2D    ! Load for 12 months on dynamics global grid
-    REAL, DIMENSION(klon_glo,12)          :: load_glo1D    ! -"- on physical global grid
+    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: psurf_glo2D
+    REAL, ALLOCATABLE, DIMENSION(:,:)     :: psurf_glo1D
+    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: load_glo2D
+    REAL, ALLOCATABLE, DIMENSION(:,:)     :: load_glo1D
     REAL, DIMENSION(iim,jjm+1)            :: vartmp
     REAL, DIMENSION(iim)                  :: lon_src              ! longitudes in file
@@ -273,4 +287,25 @@
           CALL abort_gcm('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
        END IF
+
+! 1.5) Check number of month in file opened
+!
+!**************************************************************************************************
+       ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
+!       ierr = nf90_inq_dimlen(ncid, dimid, nbr_tsteps)
+       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps) )
+       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
+         CALL abort_gcm('get_aero_fromfile', 'not the right number of months in aerosol file read',1)
+       ENDIF
+       write(lunout,*)'get_aero_fromfile: ', nbr_tsteps,' months to read from aerosols file'
+       call bcast(nbr_tsteps)
+
+! 1.6) Allocation of some variables once the number of months in the file has been determined
+!
+!**************************************************************************************************
+       ALLOCATE(psurf_glo2D(iim, jjm+1, nbr_tsteps), load_glo2D(iim, jjm+1, nbr_tsteps), stat=ierr)
+       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 0.5',1)
+       ALLOCATE(psurf_glo1D(klon_glo, nbr_tsteps), load_glo1D(klon_glo, nbr_tsteps), stat=ierr)
+       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 0.5',1)
+
 
 ! 2) Check if old or new file is avalabale.
@@ -301,5 +336,5 @@
        
      ! Allocate variables depending on the number of vertical levels
-       ALLOCATE(varmth(iim, jjm+1, klev_src), varyear(iim, jjm+1, klev_src, 12), stat=ierr)
+       ALLOCATE(varmth(iim, jjm+1, klev_src), varyear(iim, jjm+1, klev_src, nbr_tsteps), stat=ierr)
        IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 1',1)
 
@@ -323,5 +358,5 @@
           CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)) )
           
-! ++) Read surface pression, 12 month in one variable
+! ++) Read surface pression, nbr_tsteps month in one variable
 !****************************************************************************************
           ! Get variable id
@@ -363,5 +398,5 @@
 ! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
 !****************************************************************************************
-          DO imth=1, 12
+          DO imth=1, nbr_tsteps
              IF (imth.EQ.1) THEN
                 cvar=TRIM(varname)//'JAN'
@@ -424,5 +459,5 @@
           
           ! Inverse vertical levels for varyear 
-          DO imth=1, 12
+          DO imth=1, nbr_tsteps
              varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
              DO k=1, klev_src
@@ -455,5 +490,5 @@
 
 !     - Invert latitudes if necessary
-       DO imth=1, 12
+       DO imth=1, nbr_tsteps
           IF (invert_lat) THEN
 
@@ -500,5 +535,5 @@
        END DO ! imth
        
-       ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
+       ALLOCATE(varyear_glo1D(klon_glo, klev_src, nbr_tsteps), stat=ierr)
        IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 3',1)
        
@@ -532,11 +567,18 @@
     ! Allocate space for output pointer variable at local process
     IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
-    ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
-    IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5',1)
+    write(lunout,*)'get_aero_fromfile: klon, ...',klon, klev_src, nbr_tsteps 
+    ALLOCATE(pt_year(klon, klev_src, nbr_tsteps), stat=ierr)
+!    IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5.1',1)
+    IF (ASSOCIATED(pt_psurf_out)) DEALLOCATE(pt_psurf_out)
+    ALLOCATE(pt_psurf_out(klon, nbr_tsteps), stat=ierr)
+!    IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5.2',1)
+    IF (ASSOCIATED(pt_load_out)) DEALLOCATE(pt_load_out)
+    ALLOCATE(pt_load_out(klon, nbr_tsteps), stat=ierr)
+!    IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5.3',1)
 
     ! Scatter global field to local domain at local process
     CALL scatter(varyear_glo1D, pt_year)
-    CALL scatter(psurf_glo1D, psurf_out)
-    CALL scatter(load_glo1D,  load_out)
+    CALL scatter(psurf_glo1D, pt_psurf_out)
+    CALL scatter(load_glo1D,  pt_load_out)
 
 ! 7) Test for negative values
