!
MODULE readdim2
!---------------------------------------------------------------------
!- $Header: /home/ssipsl/CVSREP/sechiba/readdim2.f90,v 1.16 2000/05/11 14:48:27 ssipsl Exp $
!-
   USE ioipsl
   USE weather
!-
   IMPLICIT NONE
!-
   PRIVATE
   PUBLIC  :: forcing_read, forcing_info, forcing_grid
!-
   INTEGER, SAVE                            :: iim_full, jjm_full, llm_full, ttm_full
   INTEGER, SAVE                            :: iim_zoom, jjm_zoom, ist_zoom, jst_zoom
   REAL, SAVE, ALLOCATABLE, DIMENSION(:,:)  :: data_full, lon_full, lat_full
   REAL, SAVE, ALLOCATABLE, DIMENSION(:)    :: lev_full
   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: itau, i_index, j_index
   INTEGER, SAVE                            :: i_test, j_test
   LOGICAL, SAVE                            :: allow_weathergen, weathergen, interpol
   REAL, SAVE                               :: limit_west, limit_east, limit_north, limit_south
   REAL, SAVE                               :: merid_res, zonal_res
!-
CONTAINS
!-
!=====================================================================
!-
SUBROUTINE forcing_info(filename, iim, jjm, llm, tm, date0, dt_force,&
     & force_id)
  
!---------------------------------------------------------------------
!
!- This subroutine will get all the info from the forcing file and
!- prepare for the zoom if needed.
!- It returns to the caller the sizes of the data it will receive at
!- the forcing_read call. This is important so that the caller can
!- allocate the right space.
!-
!- filename   : name of the file to be opened
!- iim        : size in x of the forcing data
!- jjm        : size in y of the forcing data
!- llm        : number of levels in the forcing data (not yet used)
!- tm         : Time dimension of the forcing
!- date0      : The date at which the forcing file starts (julian days)
!- dt-force   : time-step of the forcing file in seconds
!- force_id   : ID of the forcing file
!-
!- ARGUMENTS
!-
  CHARACTER(LEN=*) :: filename
  INTEGER          :: iim, jjm, llm, tm
  REAL             :: date0, dt_force
  INTEGER          :: force_id
!- LOCAL
  REAL, ALLOCATABLE, DIMENSION(:,:) :: lon, lat
  REAL :: juld_1, juld_2
  LOGICAL :: debug = .FALSE.
!-
  CALL flininfo(filename,  iim_full, jjm_full, llm_full, ttm_full, force_id)
!-
  if ( debug ) WRITE(*,*) 'Details from forcing file :', iim_full, jjm_full, llm_full, ttm_full
!-
  ALLOCATE(lon(iim_full, jjm_full), lat(iim_full, jjm_full))
  ALLOCATE(itau(ttm_full))
  ALLOCATE(data_full(iim_full, jjm_full),lon_full(iim_full, jjm_full),&
       & lat_full(iim_full, jjm_full))
  ALLOCATE(lev_full(llm_full))
  ALLOCATE(i_index(iim_full), j_index(jjm_full))
!-
  CALL flinopen &
       &  (filename, .FALSE., iim_full, jjm_full, llm_full, lon_full, lat_full, &
       &   lev_full, ttm_full, itau, date0, dt_force, force_id)
!-
!- What are the alowed options for the temportal interpolation
!-
  !Config  Key  = ALLOW_WEATHERGEN
  !Config  Desc = Allow weather generator to create data
  !Config  Def  = n
  !Config  Help = This flag allows the forcing-reader to generate
  !Config         synthetic data if the data in the file is too sparse
  !Config         and the temporal resolution would not be enough to
  !Config         run the model.
!-
  allow_weathergen = .FALSE.
  CALL getin('ALLOW_WEATHERGEN',allow_weathergen)
!-
  IF (ttm_full .GE. 2) THEN
     juld_1 = itau2date(itau(1), date0, dt_force)
     juld_2 = itau2date(itau(2), date0, dt_force)
  ELSE
     WRITE(*,*) 'forcing_info : What is that only one time step in the forcing file '
     WRITE(*,*) 'That can not be right '
     STOP
  ENDIF
!- What is the distance between the two first states. From this we will deduce what is
!- to be done.
  weathergen = .FALSE.
  interpol = .FALSE.
!-
  IF ( ABS(ABS(juld_2-juld_1)-30.) .LE. 2.) THEN
    IF ( allow_weathergen ) THEN
      weathergen = .TRUE.
      WRITE(*,*) 'Using weather generator.' 
    ELSE
      WRITE(*,*) 'This seems to be a monthly file.'
      WRITE(*,*) 'We should use a weather generator with this file.'
      WRITE(*,*) 'This should be allowed in the run.def'
      STOP
    ENDIF
  ELSEIF ( ABS(juld_1-juld_2) .LE. 1./4.) THEN
    interpol = .TRUE.
    WRITE(*,*) 'We will interpolate between the forcing data time steps.' 
  ELSE
    ! Using the weather generator with data other than monthly ones probably
    ! needs some thinking.
    WRITE(*,*) 'We cannot do anything with these forcing data.'
    WRITE(*,*) 'The time step is not suitable:',ABS(juld_1-juld_2),' days.'
    STOP
  ENDIF
!-
!- redefine the forcing time step if the weather generator is activated
!-
  IF ( weathergen ) THEN
    !Config  Key  = DT_WEATHGEN
    !Config  Desc = Calling frequency of weather generator (s)
    !Config  If = ALLOW_WEATHERGEN
    !Config  Def  = 1800.
    !Config  Help = Determines how often the weather generator
    !Config         is called (time step in s). Should be equal
    !Config         to or larger than Sechiba's time step (say,
    !Config         up to 6 times Sechiba's time step or so).
    dt_force = 1800.
    CALL getin('DT_WEATHGEN',dt_force)
  ENDIF 
!-
!- Define the zoom
!-
  !Config  Key  = LIMIT_WEST
  !Config  Desc = Western limit of region
  !Config  Def  = -180.
  !Config  Help = Western limit of the region we are 
  !Config         interested in. Between -180 and +180 degrees
  !Config         The model will use the smalest regions from
  !Config         region specified here and the one of the forcing file.
!- 
  limit_west = -180.
  CALL getin('LIMIT_WEST',limit_west)
!-
  !Config  Key  = LIMIT_EAST
  !Config  Desc = Eastern limit of region
  !Config  Def  = 180.
  !Config  Help = Eastern limit of the region we are
  !Config         interested in. Between -180 and +180 degrees
  !Config         The model will use the smalest regions from
  !Config         region specified here and the one of the forcing file.
!-
  limit_east = 180.
  CALL getin('LIMIT_EAST',limit_east)
!-
  !Config  Key  = LIMIT_NORTH
  !Config  Desc = Northern limit of region
  !Config  Def  = 90.
  !Config  Help = Northern limit of the region we are
  !Config         interested in. Between +90 and -90 degrees
  !Config         The model will use the smalest regions from
  !Config         region specified here and the one of the forcing file.
!-
  limit_north = 90.
  CALL getin('LIMIT_NORTH',limit_north)
!-
  !Config  Key  = LIMIT_SOUTH
  !Config  Desc = Southern limit of region
  !Config  Def  = -90.
  !Config  Help = Southern limit of the region we are
  !Config         interested in. Between 90 and -90 degrees
  !Config         The model will use the smalest regions from
  !Config         region specified here and the one of the forcing file.
!-
  limit_south = -90.
  CALL getin('LIMIT_SOUTH',limit_south)
!- 
!- Calculate domain size
!-
  IF ( interpol ) THEN
!-
!- If we use temporal interpolation, then we cannot change the resolution (yet?)
!-
    CALL domain_size (limit_west, limit_east, limit_north, limit_south,&
         & iim_full, jjm_full, lon_full, lat_full, iim_zoom, jjm_zoom,&
         & i_index, j_index)
!-
!- If we use the weather generator, then we read zonal and meridional resolutions
!- This should be unified one day...
!-
  ELSEIF ( weathergen ) THEN
!-
    !Config  Key  = MERID_RES
    !Config  Desc = North-South Resolution
    !Config  Def  = 2.
    !Config  If = ALLOW_WEATHERGEN
    !Config  Help = North-South Resolution of the region we are
    !Config         interested in. In degrees
!-
    merid_res = 2.
    CALL getin('MERID_RES',merid_res)
!-
    !Config  Key  = ZONAL_RES
    !Config  Desc = East-West Resolution
    !Config  Def  = 2.
    !Config  If = ALLOW_WEATHERGEN
    !Config  Help = East-West Resolution of the region we are
    !Config         interested in. In degrees
!-
    zonal_res = 2.
    CALL getin('ZONAL_RES',zonal_res)
!-
    CALL weathgen_domain_size (limit_west,limit_east,limit_north,limit_south, &
                               zonal_res,merid_res,iim_zoom,jjm_zoom)
!-
!- Number of time steps is meaningless in this case
!-
    ttm_full = HUGE( ttm_full )
!-
  ELSE
!-
    STOP 'ERROR: Neither interpolation nor weather generator is specified.'
!-
  ENDIF
!-
!- Transfer the right information to the caller
!-
  iim = iim_zoom
  jjm = jjm_zoom
  llm = 1
  tm = ttm_full
!-
END SUBROUTINE forcing_info
!-
!-
!=====================================================================
SUBROUTINE forcing_read &
  & (filename, rest_id, lrstread, lrstwrite, &
  &  itauin, istp, itau_split, split, nb_spread, netrad_cons, date0,   &
  &  dt_force, iim, jjm, lon, lat, ttm, &
  &  swdown, precip, snowf, tair, &
  &  u, v, qair, pb, lwdown, &
  &  kindex, nbindex, force_id)
!---------------------------------------------------------------------
!- filename   : name of the file to be opened
!- rest_id    : ID of restart file
!- lrstread   : read restart file?
!- lrstwrite  : write restart file?
!- itauin     : time step for which we need the data
!- istp       : time step for restart file
!- itau_split : Where are we within the splitting
!-              of the time-steps of the forcing files
!-              (it decides IF we READ or not)
!- split      : The number of time steps we do
!-              between two time-steps of the forcing
!- nb_spread  : Over how many time steps do we spread the precipitation
!- netrad_cons: flag that decides if net radiation should be conserved.
!- date0      : The date at which the forcing file starts (julian days)
!- dt_force   : time-step of the forcing file in seconds
!- iim        : Size of the grid in x
!- jjm        : size of the grid in y
!- lon        : Longitudes
!- lat        : Latitudes
!- ttm        : number of time steps in all in the forcing file
!- swdown     : Downward solar radiation (W/m^2)
!- precip     : Precipitation (Rainfall) (kg/m^2s)
!- snowf      : Snowfall (kg/m^2s)
!- tair       : 2m air temperature (K)
!- u and v    : 2m (in theory !) wind speed (m/s)
!- qair       : 2m humidity (kg/kg)
!- pb         : Surface pressure (Pa)
!- lwdown     : Downward long wave radiation (W/m^2)
!- kindex     : Index of all land-points in the data
!-              (used for the gathering)
!- nbindex    : Number of land points
!- force_id   : FLINCOM file id.
!-              It is used to close the file at the end of the run.
!---------------------------------------------------------------------
   IMPLICIT NONE
!-
   CHARACTER(LEN=*) :: filename
   INTEGER :: rest_id
   LOGICAL :: lrstread, lrstwrite
   INTEGER :: itauin, istp, itau_split, split, nb_spread
   LOGICAL :: netrad_cons
   REAL    :: date0, dt_force
   INTEGER :: iim, jjm, ttm
   REAL,DIMENSION(iim,jjm) :: &
  &  lon, lat, swdown, precip, snowf, tair, u, v, qair, pb, lwdown
   INTEGER,DIMENSION(iim*jjm) :: kindex
   INTEGER :: nbindex, force_id
!-
   INTEGER :: i,j
   REAL, SAVE :: un_an, un_jour = 0.0
!-
!- Initialize if necessary
!
   IF ( un_jour .LE. 0 ) THEN
      CALL ioget_calendar(un_an, un_jour)
   ENDIF
!
   IF ( interpol ) THEN
!-
     CALL forcing_read_interpol &
         (filename, itauin, itau_split, split, nb_spread, netrad_cons, date0,   &
          dt_force, iim, jjm, lon, lat, ttm, swdown, precip, snowf, tair, &
          u, v, qair, pb, lwdown, kindex, nbindex, force_id)
!-
   ELSEIF ( weathergen ) THEN
!-
     CALL weathgen_main (itauin, istp, filename, force_id, iim, jjm, 1, &
                         rest_id, lrstread, lrstwrite, &
                         limit_west, limit_east, limit_north, limit_south, &
                         zonal_res, merid_res, lon, lat, date0, &
                         kindex, nbindex, dt_force, &
                         swdown, precip, snowf, tair, u, v, qair, pb, lwdown)
!-
!-  Transform from our old convention to the ALMA one.
!-  (Temporary fix)
!-
     DO i=1,iim
        DO j=1,jjm
           precip(i,j) = precip(i,j)/un_jour
           snowf(i,j) = snowf(i,j)/un_jour
           tair(i,j) = tair(i,j)-273.15
           pb(i,j) = pb(i,j)*100.
        ENDDO
     ENDDO
!-
   ELSE
!-
     STOP 'ERROR: Neither interpolation nor weather generator is specified.'
!-
   ENDIF
!
!-------------------------
END SUBROUTINE forcing_read
!=====================================================================
!-
!-
!=====================================================================
SUBROUTINE forcing_read_interpol &
  & (filename, itauin, itau_split, split, nb_spread, netrad_cons, date0,   &
  &  dt_force, iim, jjm, lon, lat, ttm, swdown, precip, snowf, tair, &
  &  u, v, qair, pb, lwdown, kindex, nbindex, force_id)
!---------------------------------------------------------------------
!- filename   : name of the file to be opened
!- itauin     : time step for which we need the data
!- itau_split : Where are we within the splitting
!-              of the time-steps of the forcing files
!-              (it decides IF we READ or not)
!- split      : The number of time steps we do
!-              between two time-steps of the forcing
!- nb_spread  : Over how many time steps do we spread the precipitation
!- netrad_cons: flag that decides if net radiation should be conserved.
!- date0      : The date at which the forcing file starts (julian days)
!- dt_force   : time-step of the forcing file in seconds
!- iim        : Size of the grid in x
!- jjm        : size of the grid in y
!- lon        : Longitudes
!- lat        : Latitudes
!- ttm        : number of time steps in all in the forcing file
!- swdown     : Downward solar radiation (W/m^2)
!- precip     : Precipitation (kg/m^2s)
!- snowf      : Snow fall (kg/m^2s)
!- tair       : 2m air temperature (K)
!- u and v    : 2m (in theory !) wind speed (m/s)
!- qair       : 2m humidity (kg/kg)
!- pb         : Surface pressure (Pa)
!- lwdown     : Downward long wave radiation (W/m^2)
!- kindex     : Index of all land-points in the data
!-              (used for the gathering)
!- nbindex    : Number of land points
!- force_id   : FLINCOM file id.
!-              It is used to close the file at the end of the run.
!---------------------------------------------------------------------
   IMPLICIT NONE
!-
   INTEGER,PARAMETER :: lm=1
!-
   CHARACTER(LEN=*) :: filename
   INTEGER :: itauin, itau_split, split, nb_spread
   LOGICAL :: netrad_cons
   REAL    :: date0, dt_force
   INTEGER :: iim, jjm, ttm
   REAL,DIMENSION(iim,jjm) :: &
  &  lon, lat, swdown, precip, snowf, tair, u, v, qair, pb, lwdown
   INTEGER,DIMENSION(iim*jjm) :: kindex
   INTEGER :: nbindex, force_id
!-
!!   INTEGER :: tiim, tjjm, tllm
   INTEGER, SAVE :: last_read=0, last_split
   INTEGER, SAVE :: itau_read, itau_read_n=0, itau_read_np1=0
   REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: &
  &  swdown_n, precip_n, snowf_n, tair_n, u_n, v_n, qair_n, pb_n, lwdown_n, &
  &  swdown_np1, precip_np1, snowf_np1, tair_np1, u_np1, v_np1, qair_np1,     &
  &  pb_np1, lwdown_np1, mean_sinang, sinang
   REAL, SAVE :: julian_for, un_jour, un_an
   REAL :: julian, jur, ss, rw
   INTEGER :: yy, mm, dd, is
   LOGICAL, SAVE :: initialized = .FALSE.
   LOGICAL :: check=.FALSE.
!---------------------------------------------------------------------
   IF (check) THEN
     WRITE(*,*) &
    &  " ISLSCP READ : itau_read, itau_split : ",itauin,itau_split
   ENDIF
!-
   itau_read = itauin
!-
!- This part initializes the reading of the forcing. As you can see
!- we only go through here if both time steps are zero.
!- 
   IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN
!-
     IF (check) WRITE(*,*) 'ALLOCATE all the memory needed'
!-
     ALLOCATE &
    &  (swdown_n(iim,jjm), precip_n(iim,jjm), snowf_n(iim,jjm), &
    &   tair_n(iim,jjm), u_n(iim,jjm), v_n(iim,jjm), qair_n(iim,jjm), &
    &   pb_n(iim,jjm), lwdown_n(iim,jjm))
     ALLOCATE &
    &  (swdown_np1(iim,jjm), precip_np1(iim,jjm), snowf_np1(iim,jjm), &
    &   tair_np1(iim,jjm), u_np1(iim,jjm), v_np1(iim,jjm), qair_np1(iim,jjm), &
    &   pb_np1(iim,jjm), lwdown_np1(iim,jjm))
     ALLOCATE &
    &  (mean_sinang(iim,jjm), sinang(iim,jjm))
!-
     IF (check) WRITE(*,*) 'Memory ALLOCATED'
!-
!- We need for the driver surface air temperature and humidity before the
!- the loop starts. Thus we read it here.
!-
     CALL flinget (force_id,'Tair',iim_full, jjm_full, 1, ttm,  1, 1, data_full)
     CALL forcing_zoom(data_full, tair)
     CALL flinget (force_id,'PSurf',iim_full, jjm_full, 0, ttm,  1, 1, data_full)
     CALL forcing_zoom(data_full, pb)
     CALL flinget (force_id,'Qair',iim_full, jjm_full, 1, ttm,  1, 1, data_full)
     CALL forcing_zoom(data_full, qair)
!---
     IF (check) WRITE(*,*) 'tair, psol and qair have been extracted'
!---
!--- Create the index table
!---
     CALL forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test)

     CALL ioget_calendar(un_an, un_jour)
!---
   ENDIF
!---
   IF (check) THEN
      WRITE(*,*) &
           & 'The dates : ',itau_read,itau_split,itau_read_n,itau_read_np1
   ENDIF
!---
!--- Here we do the work in case only interpolation is needed.
!---
   IF ( initialized .AND. interpol ) THEN
!---
      IF (itau_read+1 /= last_read) THEN
!---
         IF (itau_read_n == 0) THEN
            CALL flinget (force_id,'SWdown', iim_full, jjm_full, 0,ttm, &
                 &              itau_read,itau_read, data_full)
            CALL forcing_zoom(data_full, swdown_np1)
            CALL flinget (force_id,'Rainf',    iim_full, jjm_full, 0,ttm, &
                 &              itau_read,itau_read, data_full)
            CALL forcing_zoom(data_full, precip_np1)
            CALL flinget (force_id,'Snowf',    iim_full, jjm_full, 0,ttm, &
                 &              itau_read,itau_read, data_full)
            CALL forcing_zoom(data_full, snowf_np1)
            CALL flinget (force_id,'Tair',    iim_full, jjm_full, 1,ttm, &
                 &              itau_read,itau_read, data_full)
            CALL forcing_zoom(data_full, tair_np1)
            CALL flinget (force_id,'Wind_N',       iim_full, jjm_full, 1,ttm, &
                 &              itau_read,itau_read, data_full)
            CALL forcing_zoom(data_full, u_np1)
            CALL flinget (force_id,'Wind_E',       iim_full, jjm_full, 1,ttm, &
                 &              itau_read,itau_read, data_full)
            CALL forcing_zoom(data_full, v_np1)
            CALL flinget (force_id,'Qair',    iim_full, jjm_full, 1,ttm, &
                 &              itau_read,itau_read, data_full)
            CALL forcing_zoom(data_full, qair_np1)
            CALL flinget (force_id,'PSurf',    iim_full, jjm_full, 0,ttm, &
                 &              itau_read,itau_read, data_full)
            CALL forcing_zoom(data_full, pb_np1)
            CALL flinget (force_id,'LWdown',   iim_full, jjm_full, 0,ttm, &
                 &              itau_read,itau_read, data_full)
            CALL forcing_zoom(data_full, lwdown_np1)
            itau_read_n = itau_read
            itau_read_np1 = itau_read
         ENDIF
!-----
         swdown_n(:,:) = swdown_np1(:,:)
         precip_n(:,:) = precip_np1(:,:)
         tair_n(:,:)   = tair_np1(:,:)
         u_n(:,:)      = u_np1(:,:)
         v_n(:,:)      = v_np1(:,:)
         qair_n(:,:)   = qair_np1(:,:)
         pb_n(:,:)     = pb_np1(:,:)
         lwdown_n(:,:) = lwdown_np1(:,:)
         itau_read_n = itau_read_np1
!-----
         itau_read_np1 = itau_read+1
         IF (itau_read_np1 > ttm) THEN
            WRITE(*,*) 'WARNING --WARNING --WARNING --WARNING '
            WRITE(*,*) &
                 &  'WARNING : We are going back to the start of the file'
            itau_read_np1 =1
         ENDIF
!-----
!----- Get a reduced julian day !
!----- This is needed because we lack the precision on 32 bit machines.
!-----
         julian_for = itau2date(itau_read-1, date0, dt_force)
         CALL ju2ymds (julian_for, yy, mm, dd, ss)
         CALL ymds2ju (yy,1,1,0.0, jur)
         julian_for = julian_for-jur
!-----
         IF (check) THEN
            WRITE(*,*) 'Forcing for Julian day ',julian_for,'is read'
         ENDIF
!-----
         CALL flinget (force_id,'SWdown', iim_full, jjm_full, 0,ttm, &
              &              itau_read_np1,itau_read_np1, data_full)
         CALL forcing_zoom(data_full, swdown_np1)
         CALL flinget (force_id,'Rainf',    iim_full, jjm_full, 0,ttm, &
              &              itau_read_np1,itau_read_np1, data_full)
         CALL forcing_zoom(data_full, precip_np1)
         CALL flinget (force_id,'Snowf',    iim_full, jjm_full, 0,ttm, &
              &              itau_read_np1,itau_read_np1, data_full)
         CALL forcing_zoom(data_full, snowf_np1)
         CALL flinget (force_id,'Tair',    iim_full, jjm_full, 1,ttm, &
              &              itau_read_np1,itau_read_np1, data_full)
         CALL forcing_zoom(data_full, tair_np1)
         CALL flinget (force_id,'Wind_N',       iim_full, jjm_full, 1,ttm, &
              &              itau_read_np1,itau_read_np1, data_full)
         CALL forcing_zoom(data_full, u_np1)
         CALL flinget (force_id,'Wind_E',       iim_full, jjm_full, 1,ttm, &
              &              itau_read_np1,itau_read_np1, data_full)
         CALL forcing_zoom(data_full, v_np1)
         CALL flinget (force_id,'Qair',    iim_full, jjm_full, 1,ttm, &
              &              itau_read_np1,itau_read_np1, data_full)
         CALL forcing_zoom(data_full, qair_np1)
         CALL flinget (force_id,'PSurf',    iim_full, jjm_full, 0,ttm, &
              &              itau_read_np1,itau_read_np1, data_full)
         CALL forcing_zoom(data_full, pb_np1)
         CALL flinget (force_id,'LWdown',   iim_full, jjm_full, 0,ttm, &
              &              itau_read_np1,itau_read_np1, data_full)
         CALL forcing_zoom(data_full, lwdown_np1)
!-----
         last_read = itau_read_np1
!-----
!----- Compute mean solar angle for the comming period
!-----
         IF (check) WRITE(*,*) 'Going into  solarang', split, un_jour
!-----
         mean_sinang(:,:) = 0.0
         DO is=1,split
            julian = julian_for+(FLOAT(is-1)/split)*dt_force/un_jour
            CALL solarang (julian, iim, jjm, lon, lat, sinang)
            mean_sinang(:,:) = mean_sinang(:,:)+sinang(:,:)
         ENDDO
         mean_sinang(:,:) = mean_sinang(:,:)/split
!-----
      ENDIF
!---
!--- Do the interpolation
      IF (check) WRITE(*,*) 'Doing the interpolation between time steps'
!---
      rw = REAL(itau_split-1)
      qair(:,:) = ((qair_np1(:,:)-qair_n(:,:))/split)*rw+qair_n(:,:)
      tair(:,:) = ((tair_np1(:,:)-tair_n(:,:))/split)*rw+tair_n(:,:)
      pb(:,:) = ((pb_np1(:,:)-pb_n(:,:))/split)*rw+pb_n(:,:)
      u(:,:)  = ((u_np1(:,:)-u_n(:,:))/split)*rw+u_n(:,:)
      v(:,:)  = ((v_np1(:,:)-v_n(:,:))/split)*rw+v_n(:,:)
!---
!--- Here we need to allow for an option
!--- where radiative energy is conserved
!---
      IF ( netrad_cons ) THEN
         lwdown(:,:) = lwdown_n(:,:)
      ELSE
         lwdown(:,:) = ((lwdown_np1(:,:)-lwdown_n(:,:))/split)*rw &
              &             +lwdown_n(:,:)
      ENDIF
!---
!--- For the solar radiation we decompose the mean value
!--- using the zenith angle of the sun if the time step in the forcing data is 
!---- more than an hour. Else we use the standard linera interpolation
!----
      IF (check) WRITE(*,*) 'Ready to deal with the solar radiation'
!----
      IF ( dt_force .GT. 3600. ) THEN
!---
         IF ( netrad_cons ) THEN
            WRITE(*,*) 'Solar radiation can not be conserved with a timestep of ', dt_force
         ENDIF
!---
         julian = julian_for+(rw/split)*dt_force/un_jour
         IF (check) THEN
            WRITE(*,'(a,f20.10,2I3)') &
                 &  'JULIAN BEFORE SOLARANG : ',julian,itau_split,split
         ENDIF
!---
         CALL solarang(julian, iim, jjm, lon, lat, sinang)
!---
         WHERE (mean_sinang > 0.)
            swdown = swdown_n *sinang/mean_sinang
         ELSEWHERE
            swdown = 0.0
         END WHERE
!---
         WHERE (swdown > 2000. )
            swdown = 2000.
         END WHERE
      ELSE
!---
         IF ( netrad_cons ) THEN
            swdown(:,:) = swdown_n(:,:)
         ELSE
            swdown(:,:) = ((swdown_np1(:,:)-swdown_n(:,:))/split)*rw &
                 &             +swdown_n(:,:)
         ENDIF
!---
      ENDIF
      IF (check) THEN
         WRITE(*,*) 'SWNET at ',itau_split,' : ',swdown_n(i_test, j_test), &
              &           ' < ', swdown(i_test, j_test), ' < ', swdown_np1(i_test, j_test)
         WRITE(*,*) 'TAIR at  ',itau_split,' :',tair_n(i_test, j_test), &
              &           ' < ', tair(i_test, j_test), ' < ', tair_np1(i_test, j_test)
         WRITE(*,*) 'QAIR at  ',itau_split,' :',qair_n(i_test, j_test), &
              &           ' < ', qair(i_test, j_test), ' < ', qair_np1(i_test, j_test)
      ENDIF
!---
!--- For precip we suppose that the rain
!--- is the summ over the next 6 hours
!---
      IF (itau_split <= nb_spread) THEN
         precip(:,:) = precip_n(:,:)/nb_spread
      ELSE
         precip(:,:) = 0.0
      ENDIF
!---
   ENDIF
!---
!--- Here we might put the call to the weather generator ... one day.
!--- Pour le moment, le branchement entre interpolation et generateur de temps 
!--- est fait au-dessus.
!---
!-   IF ( initialized .AND. weathergen ) THEN
!-      ....
!-   ENDIF
!---
!--- At this point the code should be initialized. If not we have a problem !
!---
   IF ( (itau_read == 0).AND.(itau_split == 0) ) THEN
!---
      initialized = .TRUE.
!---
   ELSE
      IF ( .NOT. initialized ) THEN
         WRITE(*,*) 'Why is the code forcing_read not initialized ?'
         WRITE(*,*) 'Have you called it with both time-steps set to zero ?'
         STOP
      ENDIF
   ENDIF
!
!-------------------------
END SUBROUTINE forcing_read_interpol
!=====================================================================
!-
SUBROUTINE forcing_landind(iim, jjm, tair, check, nbindex, kindex, i_test, j_test)
!---
!--- This subroutine finds the indices of the land points over which the land
!--- surface scheme is going to run.
!---
!- ARGUMENTS
  INTEGER :: iim, jjm, nbindex, i_test, j_test
  LOGICAL :: check
  INTEGER :: kindex(nbindex)
  REAL    :: tair(iim,jjm)
!-
!- LOCAL
  INTEGER :: i, j, ig
!-
!-
  ig = 0
  i_test = 0
  j_test = 0
!---
  IF (MINVAL(tair(:,:)) < 100.) THEN
!----- In this case the 2m temperature is in Celsius
     DO j=1,jjm
        DO i=1,iim
           IF (tair(i,j) < 100.) THEN
              ig = ig+1
              kindex(ig) = (j-1)*iim+i
              !
              !  Here we find at random a land-point on which we can do
              !  some printouts for test.
              !
              IF (ig .GT. (iim*jjm)/2 .AND. i_test .LT. 1) THEN
                 i_test = i
                 j_test = j
                 IF (check) THEN
                    WRITE(*,*) 'The test point chosen for output is : ', i_test, j_test
                 ENDIF
              ENDIF
           ENDIF
        ENDDO
     ENDDO
  ELSE 
!----- 2m temerature is in Kelvin
     DO j=1,jjm
        DO i=1,iim
           IF (tair(i,j) < 500.) THEN
              ig = ig+1
              kindex(ig) = (j-1)*iim+i
           ENDIF
        ENDDO
     ENDDO
  ENDIF
!---
     nbindex = ig
!---
END SUBROUTINE forcing_landind
!-
!=====================================================================
!-
SUBROUTINE forcing_grid(iim,jjm,llm,lon,lat,lev)
!-
!- This subroutine calculates the longitudes and latitudes of the model grid.
!-
  INTEGER, INTENT(in)                   :: iim, jjm, llm
  REAL, DIMENSION(iim,jjm), INTENT(out) :: lon, lat
  REAL, DIMENSION(llm), INTENT(out)     :: lev
!-
  INTEGER :: i,j
  REAL :: zlev
!-
!- Should be unified one day
!-
  IF ( weathergen ) THEN
!-
    DO i = 1, iim
      lon(i,:) = limit_west + zonal_res/2. + &
                 FLOAT(i-1)*(limit_east-limit_west)/FLOAT(iim)
    ENDDO
!-
    DO j = 1, jjm
      lat(:,j) = limit_north - merid_res/2. - &
                 FLOAT(j-1)*(limit_north-limit_south)/FLOAT(jjm)
    ENDDO
!-
   !Config  Key  = HEIGHT_LEV1
   !Config  If = ALLOW_WEATHERGEN
   !Config  Desc = Height at which the atmospheric variable are given
   !Config  Def  = 2.0
   !Config  Help = The atmospheric variables (temperature and specific
   !Config         humidity are measured at a specific level.
   !Config         The height of this level is needed to compute
   !Config         corectly the turbulent transfer coefficients.
   !Config         Look at the description of the forcing
   !Config         DATA for the correct value.
!-
   zlev = 2.0
   CALL getin('HEIGHT_LEV1', zlev)
   lev(:) = zlev
!-
  ELSEIF ( interpol ) THEN
!-
    CALL forcing_zoom(lon_full, lon)
    CALL forcing_zoom(lat_full, lat)
    lev(:) = lev_full(:)
!-
  ELSE
!-
    STOP 'Neither weather generator nor temporal interpolation is specified.'
!-
  ENDIF
!-
END SUBROUTINE forcing_grid
!-
!=====================================================================
!-
SUBROUTINE forcing_zoom(x_f, x_z)
!-
!- This subroutine takes the slab of data over which we wish to run the model.
!-
REAL :: x_f(iim_full, jjm_full), x_z(iim_zoom, jjm_zoom)
!-
INTEGER :: i,j
!-
DO i=1,iim_zoom
   DO j=1,jjm_zoom
      x_z(i,j) = x_f(i_index(i),j_index(j))
   ENDDO
ENDDO
!-
END SUBROUTINE forcing_zoom
!-
!=====================================================================
SUBROUTINE solarang (julian, iim, jjm, lon, lat, csang)
!---------------------------------------------------------------------
!- This subroutine computes the solar angle according to the method
!- used by GSWP and developed by J.C. Morill.
!- See for further details :
!- http://www.atmo.arizona.edu/~morrill/swrad.html
!---------------------------------------------------------------------
   USE calendar
!-
   INTEGER,INTENT(in) :: iim, jjm
   REAL,INTENT(in)  :: julian
   REAL,DIMENSION(iim,jjm), INTENT(in)  :: lon, lat
   REAL,DIMENSION(iim,jjm), INTENT(out) :: csang
!-
   REAL :: un_an, un_jour, pi, gamma, dec, decd
   REAL :: et, gmt, le, ls, lcorr, latime, omega, omegad
   REAL :: llatd, llat
   INTEGER :: igmt, ilon, ilat
   INTEGER,SAVE,ALLOCATABLE :: zone(:), lhour(:)
!
   LOGICAL :: check = .FALSE.
!---------------------------------------------------------------------
!
   pi = 4.*ATAN(1.)
   IF (check) WRITE(*,*) 'We get the right calendar information'
   CALL ioget_calendar(un_an, un_jour)
!-
!- 1) Day angle gamma
!-
   gamma = 2.*pi*MOD(julian,un_an)/un_an
!-
!- 2) Solar declination (assumed constant for a 24 hour period)  page 7
!-    in radians
!-
   IF (check) WRITE(*,*) 'Solar declination'
    !
   dec = ( 0.006918-0.399912*COS(gamma)+0.070257*SIN(gamma) &
  &       -0.006758*COS(2*gamma)+0.000907*SIN(2*gamma)      &
  &       -0.002697*COS(3*gamma)+0.00148*SIN(3*gamma))
   decd = dec*(180/pi)
!-
!- maximum error 0.0006 rad (<3'), leads to error
!- of less than 1/2 degree in zenith angle
!-
   IF (check) WRITE(*,*) 'Equation of time'
!- 3)  Equation of time  page 11
!-
   et = ( 0.000075+0.001868*COS(gamma)-0.032077*SIN(gamma)&
  &      -0.014615*COS(2*gamma)-0.04089*SIN(2*gamma))*229.18
!-
!- Get the time zones for the current time
!-
   gmt = 24.*(julian-INT(julian))
   IF (.NOT.ALLOCATED(zone))  ALLOCATE(zone(iim))
   IF (.NOT.ALLOCATED(lhour)) ALLOCATE(lhour(iim))
!-
   igmt = NINT(gmt)
   IF (check) WRITE(*,*) 'Get time zone'
   CALL time_zone(igmt, iim, jjm, lon, zone, lhour)
!-
!- Start a loop through the grid
!-
   IF (check) WRITE(*,*) 'Start a loop through the grid'
   DO ilon=1,iim
!---
!--- 4) Local apparent time  page 13
!---
!--- ls     standard longitude (nearest 15 degree meridian)
!--- le     local longtitude
!--- lhour  local standard time
!--- latime local apparent time
!--- lcorr  longitudunal correction (minutes)
!---
     le = lon(ilon,1)
     ls = ((zone(ilon)-1)*15)-180.
     lcorr = 4.*(ls-le)*(-1)
     latime = lhour(ilon)+lcorr/60.+et/60.
     IF (latime <  0.) latime = latime+24
     IF (latime > 24.) latime = latime-24
!---
!--- 5) Hour angle omega  page 15
!---
!--- hour angle is zero at noon, postive in the morning
!--- It ranges from 180 to -180
!---
!--- omegad is omega in degrees, omega is in radians
!---
     omegad = (latime-12.)*(-15.)
     omega  = omegad*pi/180.
!---
     DO ilat=1,jjm
!-----
!----- 6)  Zenith angle  page 15
!-----
!----- csang cosine of zenith angle (radians)
!----- llatd =  local latitude in degrees
!----- llat  =  local latitude in radians
!-----
       llatd = lat(1,ilat)
       llat  = llatd*pi/180.
       csang(ilon,ilat) = &
      &  MAX(0.,SIN(dec)*SIN(llat)+COS(dec)*COS(llat)*COS(omega))
     ENDDO
   ENDDO
!----------------------
END SUBROUTINE solarang
!=====================================================================
SUBROUTINE time_zone (gmt, iim, jjm, lon, zone, lhour)
!---------------------------------------------------------------------
   INTEGER :: gmt, iim, jjm, zone(iim), lhour(iim)
   REAL :: lon(iim,jjm)
!-
   INTEGER :: deg, i, ilon, change
!---------------------------------------------------------------------
   DO ilon=1,iim
!---
!--- Convert longitude index to longtitude in degrees
!---
     deg = lon(ilon,1)
!---
!--- Determine into which time zone (15 degree interval) the
!--- longitude falls.
!---
     DO i=1,25
       IF (deg < (-187.5+(15*i))) THEN
         zone(ilon) = i
         IF (zone(ilon) == 25)   zone(ilon) = 1
         EXIT
       ENDIF
     ENDDO
!---
!--- Calculate change (in number of hours) from GMT time to
!--- local hour.  Change will be negative for zones < 13 and
!--- positive for zones > 13.
!---
!--- There is also a correction for lhour < 0 and lhour > 23
!--- to lhour between 0 and 23.
!---
     IF (zone(ilon) < 13) THEN
       change = zone(ilon)-13
       lhour(ilon) = gmt+change
     ENDIF
!---
     IF (zone(ilon) == 13) THEN
       lhour(ilon) = gmt
     ENDIF
!---
     IF (zone(ilon) > 13) THEN
       change = zone(ilon)-13
       lhour(ilon) = gmt+change
     ENDIF
     IF (lhour(ilon) <  0) lhour(ilon) = lhour(ilon)+24
     IF (lhour(ilon) > 23) lhour(ilon) = lhour(ilon)-24
!---
   ENDDO
!-----------------------
END SUBROUTINE time_zone
!=====================================================================
!
! ---------------------------------------------------------------------
!

SUBROUTINE domain_size (limit_west, limit_east, limit_north, limit_south, &
     &  iim_f, jjm_f, lon, lat, iim, jjm, iind, jind)
  
  IMPLICIT NONE
  !
  ! ARGUMENTS
  !
  REAL, INTENT(inout)  :: limit_west,limit_east,limit_north,limit_south
  INTEGER, INTENT(in)  :: iim_f, jjm_f
  REAL, INTENT(in)     :: lon(iim_f, jjm_f), lat(iim_f, jjm_f)
  INTEGER, INTENT(out) :: iim,jjm
  INTEGER, INTENT(out) :: iind(iim_f), jind(jjm_f)
  !
  ! LOCAL
  !
  INTEGER :: i, j
  REAL :: lolo
  LOGICAL :: over_dateline = .FALSE.
  !
  !
  IF ( ( ABS(limit_east) .GT. 180. ) .OR. &
       ( ABS(limit_west) .GT. 180. ) ) THEN
     WRITE(*,*) 'PROBLEME longitudes.'
     WRITE(*,*) 'Limites Ouest, Est: ',limit_west,limit_east
     STOP
  ENDIF
  !
  IF ( limit_west .GT. limit_east ) over_dateline = .TRUE.
  !
  IF ( ( limit_south .LT. -90. ) .OR. &
       ( limit_north .GT. 90. ) .OR. &
       ( limit_south .GE. limit_north ) ) THEN
     WRITE(*,*) 'PROBLEME latitudes.'
     WRITE(*,*) 'Limites Nord, Sud: ',limit_north,limit_south
     STOP
  ENDIF
  !
  ! Here we assume that the grid of the forcing data is regular. Else we would have
  ! to do more work to find the index table.
  !
  iim = 0
  DO i=1,iim_f
     !
     lolo =  lon(i,1)
     IF ( lon(i,1) .GT. 180. ) lolo =  lon(i,1) - 360.
     IF ( lon(i,1) .LT. -180. ) lolo =  lon(i,1) + 360.
     !
     IF ( over_dateline ) THEN
        IF ( lolo .LT. limit_west .OR. lolo .GT. limit_east ) THEN
           iim = iim + 1
           iind(iim) = i
        ENDIF
     ELSE
        IF ( lolo .GT. limit_west .AND. lolo .LT. limit_east ) THEN
           iim = iim + 1
           iind(iim) = i
        ENDIF
     ENDIF
     !
  ENDDO
  !
  jjm = 0
  DO j=1,jjm_f
     IF ( lat(1,j) .GT. limit_south .AND. lat(1,j) .LT. limit_north) THEN
        jjm = jjm + 1
        jind(jjm) = j
     ENDIF
  ENDDO
  !
  WRITE(*,*) 'Domain size: iim, jjm = ', iim, jjm
  
  END SUBROUTINE domain_size

!------------------
END MODULE readdim2
