Changeset 1265 for LMDZ4


Ignore:
Timestamp:
Nov 20, 2009, 5:05:39 PM (15 years ago)
Author:
Laurent Fairhead
Message:

Modifications pour utiliser le calendrier réaliste et les nouveaux fichiers a 14
valeurs mensuelles

Location:
LMDZ4/branches/LMDZ4-dev/libf/phylmd
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol.F90

    r1249 r1265  
    77CONTAINS
    88
    9 SUBROUTINE readaerosol(name_aero, type, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
     9SUBROUTINE readaerosol(name_aero, type, iyr_in, klev_src, pt_ap, pt_b, pt_out, pt_psurf, pt_load, nbr_tsteps)
    1010
    1111!****************************************************************************************
     
    3434  REAL, POINTER, DIMENSION(:)     :: pt_ap        ! Pointer for describing the vertical levels     
    3535  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
     36  REAL, POINTER, DIMENSION(:,:,:) :: pt_out       ! The massvar distributions, DIMENSION(klon,klev_src, nbr_tsteps)
     37  REAL, POINTER, DIMENSION(:,:)   :: pt_psurf     ! The massvar distributions, DIMENSION(klon, nbr_tsteps)
     38  REAL, POINTER, DIMENSION(:,:)   :: pt_load      ! The massvar distributions, DIMENSION(klon, nbr_tsteps)
     39  INTEGER                         :: nbr_tsteps   ! number of timesteps in read file (12 or 14)
    3940
    4041  ! 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.
     42  CHARACTER(len=4)                  :: cyear
     43  REAL, POINTER, DIMENSION(:,:,:)   :: pt_2
     44  REAL, POINTER, DIMENSION(:,:)     :: psurf2, load2
     45  REAL                              :: p0           ! Reference pressure
     46  INTEGER                           :: iyr1, iyr2, klev_src2
     47  INTEGER                           :: it, k, i
     48  LOGICAL, PARAMETER                :: lonlyone=.FALSE.
    4849
    4950!****************************************************************************************
     
    5657!****************************************************************************************
    5758     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)
     59     ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month
     60     ! pt_out has dimensions (klon, klev_src, nbr_tsteps)
     61     CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
    6162     
    6263
     
    6566!****************************************************************************************     
    6667     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)
     68     ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month
     69     ! pt_out has dimensions (klon, klev_src, nbr_tsteps)
     70     CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
     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, nbr_tsteps)
     79     CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
     80     
    7081     
    7182  ELSE IF (type == 'scenario') THEN
     
    7586        cyear='.nat'
    7687        WRITE(lunout,*) 'get_aero 1 iyr_in=', iyr_in,'   ',cyear
    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)
     88        ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month
     89        ! pt_out has dimensions (klon, klev_src, nbr_tsteps)
     90        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
    8091       
    8192     ELSE IF (iyr_in .GE. 2100) THEN
    8293        cyear='2100'
    8394        WRITE(lunout,*) 'get_aero 2 iyr_in=', iyr_in,'   ',cyear
    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)
     95        ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month
     96        ! pt_out has dimensions (klon, klev_src, nbr_tsteps)
     97        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
    8798       
    8899     ELSE
     
    102113        WRITE(cyear,'(I4)') iyr1
    103114        WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,'   ',cyear
    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)
     115        ! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month
     116        ! pt_out has dimensions (klon, klev_src, nbr_tsteps)
     117        CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, pt_psurf, pt_load, nbr_tsteps)
    107118       
    108119        ! If to read two decades:
     
    114125           
    115126           NULLIFY(pt_2)
    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)
     127           ! get_aero_fromfile returns pt_2 allocated and initialized with data for nbr_tsteps month
     128           ! pt_2 has dimensions (klon, klev_src, nbr_tsteps)
     129           CALL get_aero_fromfile(name_aero, cyear, klev_src2, pt_ap, pt_b, p0, pt_2, psurf2, load2, nbr_tsteps)
    119130           ! Test for same number of vertical levels
    120131           IF (klev_src /= klev_src2) THEN
     
    124135           
    125136           ! Linare interpolate to the actual year:
    126            DO it=1,12
     137           DO it=1,nbr_tsteps
    127138              DO k=1,klev_src
    128139                 DO i = 1, klon
     
    134145
    135146              DO i = 1, klon
    136                  psurf(i,it) = &
    137                       psurf(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
    138                       (psurf(i,it) - psurf2(i,it))
    139 
    140                  load(i,it) = &
    141                       load(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
    142                       (load(i,it) - load2(i,it))
     147                 pt_psurf(i,it) = &
     148                      pt_psurf(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
     149                      (pt_psurf(i,it) - psurf2(i,it))
     150
     151                 pt_load(i,it) = &
     152                      pt_load(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
     153                      (pt_load(i,it) - load2(i,it))
    143154              END DO
    144155           END DO
     
    146157           ! Deallocate pt_2 no more needed
    147158           DEALLOCATE(pt_2)
     159           DEALLOCATE(psurf2)
     160           DEALLOCATE(load2)
    148161           
    149162        END IF ! lonlyone
     
    159172
    160173
    161   SUBROUTINE get_aero_fromfile(varname, cyr, klev_src, pt_ap, pt_b, p0, pt_year, psurf_out, load_out)
    162 !****************************************************************************************
    163 ! Read 12 month aerosol from file and distribute to local process on physical grid.
     174  SUBROUTINE get_aero_fromfile(varname, cyr, klev_src, pt_ap, pt_b, p0, pt_year, pt_psurf_out, pt_load_out, nbr_tsteps)
     175!****************************************************************************************
     176! Read nbr_tsteps month aerosol from file and distribute to local process on physical grid.
    164177! Vertical levels, klev_src, may differ from model levels if new file format.
    165178!
     
    197210    REAL, POINTER, DIMENSION(:)           :: pt_b         ! Pointer for describing the vertical levels     
    198211    REAL                                  :: p0           ! Reference pressure value
    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
     212    REAL, POINTER, DIMENSION(:,:,:)       :: pt_year      ! Pointer-variabale from file, nbr_tsteps month, grid : klon,klev_src
     213    REAL, POINTER, DIMENSION(:,:)         :: pt_psurf_out    ! Surface pression for nbr_tsteps months
     214    REAL, POINTER, DIMENSION(:,:)         :: pt_load_out     ! Aerosol mass load in each column
     215    INTEGER                               :: nbr_tsteps ! number of month in file read
    202216
    203217! Local variables
     
    209223    REAL                  :: npole, spole
    210224    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varmth
    211     REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear       ! Global variable read from file, 12 month
    212     REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varyear_glo1D !(klon_glo, klev_src, 12)
     225    REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear       ! Global variable read from file, nbr_tsteps month
     226    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: varyear_glo1D !(klon_glo, klev_src, nbr_tsteps)
    213227    REAL, ALLOCATABLE, DIMENSION(:)       :: varktmp
    214228
    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
     229    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: psurf_glo2D
     230    REAL, ALLOCATABLE, DIMENSION(:,:)     :: psurf_glo1D
     231    REAL, ALLOCATABLE, DIMENSION(:,:,:)   :: load_glo2D
     232    REAL, ALLOCATABLE, DIMENSION(:,:)     :: load_glo1D
    219233    REAL, DIMENSION(iim,jjm+1)            :: vartmp
    220234    REAL, DIMENSION(iim)                  :: lon_src              ! longitudes in file
     
    273287          CALL abort_gcm('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
    274288       END IF
     289
     290! 1.5) Check number of month in file opened
     291!
     292!**************************************************************************************************
     293       ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
     294!       ierr = nf90_inq_dimlen(ncid, dimid, nbr_tsteps)
     295       CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps) )
     296       IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
     297         CALL abort_gcm('get_aero_fromfile', 'not the right number of months in aerosol file read',1)
     298       ENDIF
     299       write(lunout,*)'get_aero_fromfile: ', nbr_tsteps,' months to read from aerosols file'
     300       call bcast(nbr_tsteps)
     301
     302! 1.6) Allocation of some variables once the number of months in the file has been determined
     303!
     304!**************************************************************************************************
     305       ALLOCATE(psurf_glo2D(iim, jjm+1, nbr_tsteps), load_glo2D(iim, jjm+1, nbr_tsteps), stat=ierr)
     306       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 0.5',1)
     307       ALLOCATE(psurf_glo1D(klon_glo, nbr_tsteps), load_glo1D(klon_glo, nbr_tsteps), stat=ierr)
     308       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 0.5',1)
     309
    275310
    276311! 2) Check if old or new file is avalabale.
     
    301336       
    302337     ! Allocate variables depending on the number of vertical levels
    303        ALLOCATE(varmth(iim, jjm+1, klev_src), varyear(iim, jjm+1, klev_src, 12), stat=ierr)
     338       ALLOCATE(varmth(iim, jjm+1, klev_src), varyear(iim, jjm+1, klev_src, nbr_tsteps), stat=ierr)
    304339       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 1',1)
    305340
     
    323358          CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)) )
    324359         
    325 ! ++) Read surface pression, 12 month in one variable
     360! ++) Read surface pression, nbr_tsteps month in one variable
    326361!****************************************************************************************
    327362          ! Get variable id
     
    363398! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
    364399!****************************************************************************************
    365           DO imth=1, 12
     400          DO imth=1, nbr_tsteps
    366401             IF (imth.EQ.1) THEN
    367402                cvar=TRIM(varname)//'JAN'
     
    424459         
    425460          ! Inverse vertical levels for varyear
    426           DO imth=1, 12
     461          DO imth=1, nbr_tsteps
    427462             varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
    428463             DO k=1, klev_src
     
    455490
    456491!     - Invert latitudes if necessary
    457        DO imth=1, 12
     492       DO imth=1, nbr_tsteps
    458493          IF (invert_lat) THEN
    459494
     
    500535       END DO ! imth
    501536       
    502        ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
     537       ALLOCATE(varyear_glo1D(klon_glo, klev_src, nbr_tsteps), stat=ierr)
    503538       IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 3',1)
    504539       
     
    532567    ! Allocate space for output pointer variable at local process
    533568    IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
    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)
     569    write(lunout,*)'get_aero_fromfile: klon, ...',klon, klev_src, nbr_tsteps
     570    ALLOCATE(pt_year(klon, klev_src, nbr_tsteps), stat=ierr)
     571!    IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5.1',1)
     572    IF (ASSOCIATED(pt_psurf_out)) DEALLOCATE(pt_psurf_out)
     573    ALLOCATE(pt_psurf_out(klon, nbr_tsteps), stat=ierr)
     574!    IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5.2',1)
     575    IF (ASSOCIATED(pt_load_out)) DEALLOCATE(pt_load_out)
     576    ALLOCATE(pt_load_out(klon, nbr_tsteps), stat=ierr)
     577!    IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 5.3',1)
    536578
    537579    ! Scatter global field to local domain at local process
    538580    CALL scatter(varyear_glo1D, pt_year)
    539     CALL scatter(psurf_glo1D, psurf_out)
    540     CALL scatter(load_glo1D,  load_out)
     581    CALL scatter(psurf_glo1D, pt_psurf_out)
     582    CALL scatter(load_glo1D,  pt_load_out)
    541583
    542584! 7) Test for negative values
  • LMDZ4/branches/LMDZ4-dev/libf/phylmd/readaerosol_interp.F90

    r1249 r1265  
    1313! 4) Test for negative mass values
    1414
     15  USE ioipsl
    1516  USE dimphy, ONLY : klev,klon
    1617  USE mod_phys_lmdz_para, ONLY : mpi_rank 
     
    5051  INTEGER                         :: i, k, ierr
    5152  INTEGER                         :: iday, iyr, lmt_pas
    52   INTEGER                         :: im, day1, day2, im2
     53!  INTEGER                         :: im, day1, day2, im2
     54  INTEGER                         :: im, im2
     55  REAL                            :: day1, day2
    5356  INTEGER                         :: pi_klev_src ! Only for testing purpose
    5457  INTEGER, SAVE                   :: klev_src    ! Number of vertical levles in source field
     
    7679
    7780  REAL, DIMENSION(:,:,:), POINTER   :: pt_tmp      ! Pointer allocated in readaerosol
     81  REAL, DIMENSION(:,:), POINTER     :: pt_tmp_surf ! Pointer allocated in readaerosol
     82  REAL, DIMENSION(:,:), POINTER     :: pt_tmp_load ! Pointer allocated in readaerosol
     83  REAL, DIMENSION(:,:), POINTER     :: pt_tmp_pi_surf ! Pointer allocated in readaerosol
     84  REAL, DIMENSION(:,:), POINTER     :: pt_tmp_pi_load ! Pointer allocated in readaerosol
    7885  REAL, POINTER, DIMENSION(:), SAVE :: pt_ap, pt_b ! Pointer for describing the vertical levels
    7986!$OMP THREADPRIVATE(pt_ap, pt_b)
    80 
     87  INTEGER, SAVE                     :: nbr_tsteps ! number of time steps in file read
     88  REAL, DIMENSION(14), SAVE         :: month_len, month_start, month_mid
     89!$OMP THREADPRIVATE(month_len, month_start, month_mid)
     90  REAL                              :: jDay
     91
     92!$OMP THREADPRIVATE(nbr_tsteps)
    8193  LOGICAL            :: lnewday      ! Indicates if first time step at a new day
    8294  LOGICAL            :: OLDNEWDAY
     
    98110
    99111  ! Use phys_cal_mod
    100   !iday= day_cur
    101   !iyr = year_cur
    102   !im  = mth_cur
     112! iday= day_cur
     113! iyr = year_cur
     114! im  = mth_cur
    103115
    104116  iday = INT(r_day)
     
    107119  iyr  = iyr + annee_ref      ! year of the run   
    108120  im   = iday/30 +1           ! the actual month
     121  CALL ymds2ju(iyr, im, iday, 0., jDay)
    109122
    110123  IF(MOD(itap-1,NINT(86400./pdtphys)) == 0)THEN
     
    140153     IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 2',1)
    141154
    142      ALLOCATE( psurf_year(klon, 12, naero_spc), pi_psurf_year(klon, 12, naero_spc), stat=ierr)
    143      IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 3',1)
    144 
    145      ALLOCATE( load_year(klon, 12, naero_spc), pi_load_year(klon, 12, naero_spc), stat=ierr)
    146      IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 4',1)
     155!     ALLOCATE( psurf_year(klon, 12, naero_spc), pi_psurf_year(klon, 12, naero_spc), stat=ierr)
     156!     IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 3',1)
     157!
     158!     ALLOCATE( load_year(klon, 12, naero_spc), pi_load_year(klon, 12, naero_spc), stat=ierr)
     159!     IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 4',1)
    147160
    148161     lnewday=.TRUE.
     
    159172  IF ( (first .OR. iday==0) .AND. lnewday ) THEN
    160173     NULLIFY(pt_tmp)
     174     NULLIFY(pt_tmp_surf)
     175     NULLIFY(pt_tmp_pi_surf)
     176     NULLIFY(pt_tmp_load)
     177     NULLIFY(pt_tmp_pi_load)
    161178
    162179     ! Reading values corresponding to the closest year taking into count the choice of aer_type.
    163180     ! For aer_type=scenario interpolation between 2 data sets is done in readaerosol.
    164181     CALL readaerosol(name_aero(id_aero), aer_type, iyr, klev_src, pt_ap, pt_b, pt_tmp, &
    165           psurf_year(:,:,id_aero), load_year(:,:,id_aero))
     182          pt_tmp_surf,pt_tmp_load,nbr_tsteps)
    166183     IF (.NOT. ALLOCATED(var_year)) THEN
    167         ALLOCATE(var_year(klon, klev_src, 12, naero_spc), stat=ierr)
    168         IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5',1)
     184        ALLOCATE(var_year(klon, klev_src, nbr_tsteps, naero_spc), stat=ierr)
     185        IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5.1',1)
     186        ALLOCATE( psurf_year(klon, nbr_tsteps, naero_spc), load_year(klon, nbr_tsteps, naero_spc), stat=ierr)
     187        IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 5.2',1)
    169188     END IF
    170189     var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
     190     psurf_year(:,:,id_aero) = pt_tmp_surf(:,:)
     191     load_year(:,:,id_aero) = pt_tmp_load(:,:)
    171192
    172193     ! Reading values corresponding to the preindustrial concentrations.
    173194     CALL readaerosol(name_aero(id_aero), 'preind', iyr, pi_klev_src, pt_ap, pt_b, pt_tmp, &
    174           pi_psurf_year(:,:,id_aero), pi_load_year(:,:,id_aero))
     195          pt_tmp_surf,pt_tmp_load,nbr_tsteps)
    175196
    176197     ! klev_src must be the same in both files.
     
    183204
    184205     IF (.NOT. ALLOCATED(pi_var_year)) THEN
    185         ALLOCATE(pi_var_year(klon, klev_src, 12, naero_spc), stat=ierr)
    186         IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6',1)
     206        ALLOCATE(pi_var_year(klon, klev_src, nbr_tsteps, naero_spc), stat=ierr)
     207        IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6.1',1)
     208        ALLOCATE( pi_psurf_year(klon, nbr_tsteps, naero_spc), pi_load_year(klon, nbr_tsteps, naero_spc), stat=ierr)
     209        IF (ierr /= 0) CALL abort_gcm('readaerosol_interp', 'pb in allocation 6.2',1)
    187210     END IF
    188211     pi_var_year(:,:,:,id_aero) = pt_tmp(:,:,:)
     212     pi_psurf_year(:,:,id_aero) = pt_tmp_surf(:,:)
     213     pi_load_year(:,:,id_aero) = pt_tmp_load(:,:)
    189214   
    190215     IF (debug) THEN
    191         CALL writefield_phy('var_year_jan',var_year(:,:,1,id_aero),klev_src)
    192         CALL writefield_phy('var_year_dec',var_year(:,:,12,id_aero),klev_src)
     216        CALL writefield_phy('var_year_first',var_year(:,:,1,id_aero),klev_src)
     217        CALL writefield_phy('var_year_last',var_year(:,:,nbr_tsteps,id_aero),klev_src)
    193218        CALL writefield_phy('psurf_src',psurf_year(:,:,id_aero),1)
    194219        CALL writefield_phy('pi_psurf_src',pi_psurf_year(:,:,id_aero),1)
     
    199224     ! Pointer no more useful, deallocate.
    200225     DEALLOCATE(pt_tmp)
     226     DEALLOCATE(pt_tmp_surf)
     227     DEALLOCATE(pt_tmp_load)
    201228
    202229     ! Test if vertical interpolation will be needed.
     
    220247     END IF
    221248
     249!    Calendar initialisation
     250!
     251     DO i = 2, 13
     252       month_len(i) = float(ioget_mon_len(year_cur, i-1))
     253       CALL ymds2ju(year_cur, i-1, 1, 0.0, month_start(i))
     254     ENDDO
     255     month_len(1) = float(ioget_mon_len(year_cur-1, 12))
     256     CALL ymds2ju(year_cur-1, 12, 1, 0.0, month_start(1))
     257     month_len(14) = float(ioget_mon_len(year_cur+1, 1))
     258     CALL ymds2ju(year_cur+1, 1, 1, 0.0, month_start(14))
     259     month_mid(:) = month_start (:) + month_len(:)/2.
     260
     261     write(55,*)'month_len = ',month_len
     262     write(55,*)'month_start = ',month_start
     263     write(55,*)'month_mid = ',month_mid
     264       
    222265  END IF  ! IF ( (first .OR. iday==0) .AND. lnewday ) THEN
    223266 
     
    233276!
    234277!****************************************************************************************
    235      ! Find which months and days to use for time interpolation
    236      IF (iday < im*30-15) THEN
    237         ! in the first half of the month use month before and actual month
    238         im2=im-1
    239         day2 = im2*30-15
    240         day1 = im2*30+15
    241         IF (im2 <= 0) THEN
    242            ! the month is january, thus the month before december
    243            im2=12
    244         END IF
     278    ! Find which months and days to use for time interpolation
     279     IF (nbr_tsteps == 12) then
     280       IF (jDay < month_mid(im+1)) THEN
     281          im2=im-1
     282          day2 = month_mid(im2+1)
     283          day1 = month_mid(im+1)
     284          IF (im2 <= 0) THEN
     285             ! the month is january, thus the month before december
     286             im2=12
     287          END IF
     288       ELSE
     289          ! the second half of the month
     290          im2=im+1
     291          day2 = month_mid(im+1)
     292          day1 = month_mid(im2+1)
     293          IF (im2 > 12) THEN
     294             ! the month is december, the following thus january
     295             im2=1
     296          ENDIF
     297       END IF
     298     ELSE IF (nbr_tsteps == 14) then
     299       im = im + 1
     300       IF (jDay < month_mid(im)) THEN
     301          ! in the first half of the month use month before and actual month
     302          im2=im-1
     303          day2 = month_mid(im2)
     304          day1 = month_mid(im)
     305       ELSE
     306          ! the second half of the month
     307          im2=im+1
     308          day2 = month_mid(im)
     309          day1 = month_mid(im2)
     310       END IF
    245311     ELSE
    246         ! the second half of the month
    247         im2=im+1
    248         IF (im2 > 12) THEN
    249            ! the month is december, the following thus january
    250            im2=1
    251         ENDIF
    252         day2 = im*30-15
    253         day1 = im*30+15
    254      END IF
    255      
     312       CALL abort_gcm('readaerosol_interp', 'number of months undefined',1)
     313     ENDIF
     314
     315!     ! Find which months and days to use for time interpolation
     316!     IF (iday < im*30-15) THEN
     317!        ! in the first half of the month use month before and actual month
     318!        im2=im-1
     319!        day2 = im2*30-15
     320!        day1 = im2*30+15
     321!        IF (im2 <= 0) THEN
     322!           ! the month is january, thus the month before december
     323!           im2=12
     324!        END IF
     325!     ELSE
     326!        ! the second half of the month
     327!        im2=im+1
     328!        IF (im2 > 12) THEN
     329!           ! the month is december, the following thus january
     330!           im2=1
     331!        ENDIF
     332!        day2 = im*30-15
     333!        day1 = im*30+15
     334!     END IF
     335!      jDay = jDay+1
     336!      write(55,*)'iday, jDay, im, im2, day1, day2 = ',iday, jDay, im, im2, day1, day2, nbr_tsteps
     337 
    256338     ! Time interpolation, still on vertical source grid
    257339     ALLOCATE(tmp1(klon,klev_src), tmp2(klon,klev_src),stat=ierr)
     
    265347        DO i=1,klon
    266348           tmp1(i,k) = &
    267                 var_year(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
     349                var_year(i,k,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
    268350                (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero))
    269351           
    270352           tmp2(i,k) = &
    271                 pi_var_year(i,k,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
     353                pi_var_year(i,k,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
    272354                (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero))
    273355        END DO
     
    277359     DO i=1,klon
    278360        psurf_day(i) = &
    279              psurf_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
     361             psurf_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
    280362             (psurf_year(i,im2,id_aero) - psurf_year(i,im,id_aero))
    281363       
    282364        pi_psurf_day(i) = &
    283              pi_psurf_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
     365             pi_psurf_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
    284366             (pi_psurf_year(i,im2,id_aero) - pi_psurf_year(i,im,id_aero))
    285367     END DO
     
    288370     DO i=1,klon
    289371        load_src(i) = &
    290              load_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
     372             load_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
    291373             (load_year(i,im2,id_aero) - load_year(i,im,id_aero))
    292374       
    293375        pi_load_src(i) = &
    294              pi_load_year(i,im2,id_aero) - FLOAT(iday-day2)/FLOAT(day1-day2) * &
     376             pi_load_year(i,im2,id_aero) - FLOAT(jDay-day2)/FLOAT(day1-day2) * &
    295377             (pi_load_year(i,im2,id_aero) - pi_load_year(i,im,id_aero))
    296378     END DO
     
    447529              ! Test for var_day
    448530              IF (var_day(i,k,id_aero) < 0.) THEN
    449                  IF (iday-day2 < 0.) WRITE(lunout,*) 'iday-day2=',iday-day2
     531                 IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2
    450532                 IF (var_year(i,k,im2,id_aero) - var_year(i,k,im,id_aero) < 0.) THEN
    451533                    WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', &
     
    466548              ! Test for pi_var_day
    467549              IF (pi_var_day(i,k,id_aero) < 0.) THEN
    468                  IF (iday-day2 < 0.) WRITE(lunout,*) 'iday-day2=',iday-day2
     550                 IF (jDay-day2 < 0.) WRITE(lunout,*) 'jDay-day2=',jDay-day2
    469551                 IF (pi_var_year(i,k,im2,id_aero) - pi_var_year(i,k,im,id_aero) < 0.) THEN
    470552                    WRITE(lunout,*) trim(name_aero(id_aero)),'(i,k,im2)-', &
Note: See TracChangeset for help on using the changeset viewer.