source: perso/abdelouhab.djerrah/ORCHIDEE_OL/weather.f90 @ 855

! =================================================================================================================================
! MODULE       : weather
!
! CONTACT      : orchidee-help _at_ ipsl.jussieu.fr
!
! LICENCE      : IPSL (2006)
! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
!
!>\BRIEF        This module simulates weather.
!!
!!\n DESCRIPTION: None
!!
!! RECENT CHANGE(S): None
!!
!! REFERENCE(S) :
!!
!! SVN          :
!! $HeadURL: $
!! $Date: $
!! $Revision: $
!! \n
!_ ================================================================================================================================

MODULE weather
!---------------------------------------------------------------------
  USE netcdf
!-
  USE defprec
  USE ioipsl
  USE constantes
  USE parallel
  USE grid, ONLY : year,month,day,sec
!-
  IMPLICIT NONE
!-
  PRIVATE
  PUBLIC weathgen_main, weathgen_domain_size, weathgen_init, weathgen_read_file
!
! Only for root proc
  INTEGER, SAVE                              :: iim_file, jjm_file, llm_file, ttm_file !(unitless)
  INTEGER,DIMENSION(:,:),SAVE,ALLOCATABLE    :: ncorr !(unitless)
  INTEGER,DIMENSION(:,:,:),SAVE,ALLOCATABLE  :: icorr,jcorr !(unitless)
  INTEGER,SAVE                               :: i_cut, n_agg !(unitless)

  REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinwet   !! climatological wet days + anomaly @tex $(days.month^{-1})$ @endtex
  REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinprec  !! climatological precipition + anomaly @tex $(mm.day^{-1})$ @endtex
  REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xint     !! climatological temp + anomaly (C)
  REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinq     !! climatological relative humidity + anomaly (0-1, unitless)
  REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xinwind  !! climatological wind speed + anomaly @tex $(m.s^{-1})$ @endtex
  REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xincld   !! climatological cloudiness + anomaly (0-1, unitless)
  REAL,DIMENSION(:,:),SAVE,ALLOCATABLE :: xintrng  !! climatological temp range + anomaly (C)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: xintopo  !! topography (m)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: lat_land !! latitudes of land points (unitless)
!-
! daily values
!-
  REAL,SAVE :: julian_last
  INTEGER,DIMENSION(:),SAVE,ALLOCATABLE :: iwet    !! flag for wet day / dry day (0-1, unitless)
!-
! today's values (m0 means "minus 0")
!-
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: psurfm0  !! surface pressure today(Pa)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: cloudm0  !! cloud fraction today (0-1, unitless)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: tmaxm0   !! maximum daily temperature today (K)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: tminm0   !! minimum daily temperature today (K)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: qdm0     !! daily average specific humidity today (kg_h2o/kg_air)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: udm0     !! daily average wind speed today  @tex $(m.s^{-1})$ @endtex
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: precipm0 !! daily precitation today @tex $(mm.day^{-1})$ @endtex
!-
! yesterday's values (m1 means "minus 1")
!-
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: psurfm1  !! surface pressure yesterday (Pa)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: cloudm1  !! cloud fraction yesterday (0-1, unitless)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: tmaxm1   !! maximum daily temperature yesterday (K)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: tminm1   !! minimum daily temperature yesterday (K)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: qdm1     !! daily average specific humidity yesterday (kg_h2o/kg_air)
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: udm1     !! daily average wind speed  yesterday  @tex $(m.s^{-1})$ @endtex 
  REAL,DIMENSION(:),SAVE,ALLOCATABLE   :: precipm1 !! daily precitation  yesterday @tex $(mm.day^{-1})$ @endtex 
!-
! other
!-
  REAL,PARAMETER   :: undef=1.e+20  !! special value
  INTEGER,SAVE     :: ipprec        !! statistical (0) or prescribed (1) daily values

  LOGICAL,SAVE                  :: precip_exact  !! respect monthly precipitation
  INTEGER,DIMENSION(31,12),SAVE :: jour_precip


  INTEGER,PARAMETER      :: seedsize_max = 300  !! max size of random seed
  LOGICAL,SAVE           :: dump_weather
  CHARACTER(LEN=20),SAVE :: dump_weather_file
  LOGICAL,SAVE           :: gathered
  INTEGER,SAVE           :: dump_id

!
  REAL,PARAMETER         :: zero_t=273.16   !! Absolute zero
!-
! Parametres orbitaux:
!-
  REAL,SAVE :: ecc    !! Eccentricity
  REAL,SAVE :: perh   !! Longitude of perihelie
  REAL,SAVE :: xob    !! obliquity
!-
  INTEGER,PARAMETER :: nmon = 12  !! Number of months

!
  CHARACTER(LEN=3),DIMENSION(12) :: cal = &             !! Months of the year (unitless)
 &  (/ 'JAN','FEB','MAR','APR','MAY','JUN', &
 &     'JUL','AUG','SEP','OCT','NOV','DEC' /)
  INTEGER,DIMENSION(12),SAVE :: ndaypm = &              !! Number of days per month (noleap year) (unitless)
 &  (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /)
  INTEGER,SAVE :: soldownid, rainfid, snowfid, lwradid, &
 &                tairid, qairid,  psolid, uid, vid, &
 &                time_id, timestp_id
!-
  INTEGER,SAVE :: n_rtp = nf90_real4      !! Parameters for NETCDF 
!- 
  INTEGER  :: ALLOC_ERR                   !! Flag for dynamic allocations 
!-
  CHARACTER(LEN=20),SAVE :: calendar_str   !! Calendar type
!-
  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: kindex_w  !! Land points index
  INTEGER, SAVE                            :: nbindex_w  !! Land points index
!-
  INTEGER, PARAMETER :: termcol = 100     !! Plot of projection file => grid, number of column on terminal
!80
  CONTAINS
!-
!===
!-

!! ================================================================================================================================
!! SUBROUTINE   : daily
!!
!>\BRIEF          This routine generates daily weather conditions from monthly-mean
!! climatic parameters.
!!
!! DESCRIPTION  :  Specifically, this routine generates daily values of 
!!                 - daily total precipitation \n
!!                 - daily maximum temperature \n
!!                 - daily minimum temperature \n
!!                 - daily average cloud cover \n
!!                 - daily average relative humidity \n
!!                 - daily average wind speed \n
!!                 In order to generate daily weather conditions, the model uses
!!                 a series of 'weather generator' approaches,
!!                 which generate random combinations of weather conditions
!!                 based upon the climatological conditions in general. \n
!!                 This weather generator is based upon the so-called Richardson
!!                 weather generator. \n
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::psurf, ::cloud, ::tmax, ::tmin, ::qd, ::ud, ::precip 
!!
!! REFERENCE(S) :
!! -  Geng, S., F.W.T. Penning de Vries, and L. Supit, 1985:  A simple
!! method for generating rainfall data, Agricultural and Forest
!! Meteorology, 36, 363-376.
!! -  Richardson, C. W. and Wright, D. A., 1984: WGEN: A model for
!! generating daily weather variables: U. S. Department of
!! Agriculture, Agricultural Research Service. 
!! - Richardson, C., 1981: Stochastic simulation of daily
!! precipitation, temperature, and solar radiation. Water Resources
!! Research 17, 182-190. 
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE daily &
 &  (npoi, imonth, iday, cloud, tmax, tmin, precip, qd, ud, psurf)

!-
! in & out: global variables
!-

! INTEGER,INTENT(INOUT):: iwet(npoi) !! wet day / dry day flag!*
!-

  !! 0. Parameters and variables declaration

  REAL,PARAMETER :: rair = 287.     !! Molecular weight of one mole of dry air @tex $(g.mol{-1})$ @endtex
  REAL,PARAMETER :: grav = 9.81     !! Acceleration of the gravity @tex $(m.s^{-2})$ @endtex
  REAL,PARAMETER :: pi = 3.1415927  !! pi

  !! 0.1 Input variables

  INTEGER,INTENT(IN) :: npoi         !! total number of land points
  INTEGER,INTENT(IN) :: imonth, iday

  !! 0.2 Output variables

  REAL,INTENT(OUT) :: psurf(npoi)   !! surface pressure (Pa)
  REAL,INTENT(OUT) :: cloud(npoi)   !! cloud fraction (0-1, unitless)
  REAL,INTENT(OUT) :: tmax(npoi)    !! maximum daily temperature (K)
  REAL,INTENT(OUT) :: tmin(npoi)    !! minimum daily temperature (K)
  REAL,INTENT(OUT) :: qd(npoi)      !! daily average specific humidity (kg_h2o/kg_air)
  REAL,INTENT(OUT) :: ud(npoi)      !! daily average wind speed  @tex $(m.s^{-1})$ @endtex
  REAL,INTENT(OUT) :: precip(npoi)  !! daily precitation @tex $(mm.day^{-1})$ @endtex

  !! 0.4 Local variables

  REAL :: td(npoi)                                 !! daily average temperature (K)
  REAL,allocatable,save,dimension(:,:) ::  xstore  !! weather generator 'memory' matrix
  REAL :: ee(3), r(3), rr(npoi,3)
  REAL :: alpha(npoi), rndnum, pwd, pww
  REAL :: beta(npoi)
  REAL :: pwet(npoi)
  REAL :: rwork
  REAL :: omcloud, omqd, omtmax
  REAL :: cloudm, cloudw, cloudd
  REAL :: cloude(npoi), clouds(npoi)
  REAL :: tmaxd, tmaxw, tmaxm
  REAL :: tminm
  REAL :: tmins(npoi), tmaxs(npoi)
  REAL :: tmine(npoi), tmaxe(npoi)
  REAL :: qdm(npoi),qdd(npoi),qde(npoi),qdw(npoi),qdup(npoi),qdlow(npoi)
  INTEGER :: i,j,k
  REAL :: amn,b1,b2,b3,eud,rn,rn1,rn2,rn3,rn4,s1,s2,s12,x1,y, z(npoi)
  REAL :: aa(npoi),ab(npoi),tr1(npoi), tr2(npoi)
  REAL :: tdm,trngm,tdum
  REAL :: qsattd(npoi)
  INTEGER :: it1w, it2w
  REAL :: dt
  REAL :: rainpwd(npoi)
  INTEGER :: not_ok(npoi)
  INTEGER :: count_not_ok,count_not_ok_g
  LOGICAL,SAVE :: firstcall = .TRUE.
  INTEGER,save :: npoi0
  REAL :: xx,e

  REAL, DIMENSION(3,3) :: a,b              !! define autocorrelation matrices for Richardson generator
                                           !!
                                           !! note that this matrix should be based upon a statistical
                                           !! analysis of regional weather patterns
                                           !!
                                           !! for global simulations, we use 'nominal' values

  LOGICAL :: Warning_aa_ab(npoi), Warning_iwet(npoi) !! Warnings
!---------------------------------------------------------------------
!-
! initial setup for daily climate calculations
!-
  a(1,:) = (/ 0.600,  0.500,  0.005 /)
  a(2,:) = (/ 0.010,  0.250,  0.005 /)
  a(3,:) = (/ 0.020,  0.125,  0.250 /)
!-
  b(1,:) = (/ 0.500,  0.250, -0.250 /)
  b(2,:) = (/ 0.000,  0.500,  0.250 /)
  b(3,:) = (/ 0.000,  0.000,  0.500 /)
!-
  e = EXP(1.)
! GK240100
  IF (firstcall) THEN
    firstcall = .FALSE.
    ALLOC_ERR=-1
    ALLOCATE(xstore(npoi,3), STAT=ALLOC_ERR)
    IF (ALLOC_ERR/=0) THEN
      WRITE(numout,*) "ERROR IN ALLOCATION of xstore : ",ALLOC_ERR
      STOP 
    ENDIF
    xstore(:,:) = 0.
    npoi0 = npoi
  ELSE IF (npoi /= npoi0) THEN
    WRITE(numout,*) 'Domain size old, new: ',npoi0,npoi
    STOP 'WG Daily: Problem: Domain has changed since last call'
  ENDIF
!-
! define working variables
!-
  rwork = (grav/rair/0.0065)
!-
! 'omega' parameters used to calculate differences in expected
! climatic parameters on wet and dry days
!
! following logic of weather generator used in the EPIC crop model
!
! omcloud -- cloud cover
! omqd    -- humidity
! omtmax  -- maximum temperature
!-
  omcloud = 0.90    ! originally 0.90
  omqd    = 0.50    ! originally 0.50
  omtmax  = 0.75    ! originally 0.75
!-
! calculate weighting factors used in interpolating climatological
! monthly-mean input values to daily-mean values
!-
! this is a simple linear interpolation technique that takes into
! account the length of each month
!-
  IF (ipprec == 0) THEN
    IF (REAL(iday) < REAL(ndaypm(imonth)+1)/2.0) then
      it1w = imonth-1
      it2w = imonth
      dt   = (REAL(iday)-0.5)/ndaypm(imonth)+0.5
    ELSE
      it1w = imonth
      it2w = imonth+1
      dt   = (REAL(iday)-0.5)/ndaypm(imonth)-0.5
    ENDIF
    if (it1w <  1)  it1w = 12
    if (it2w > 12)  it2w = 1
  ELSE
     dt = -1.
     it1w = -1
     it2w = -1
  ENDIF
!-
  IF (ipprec == 0) THEN
!---
!-- use weather generator to create daily statistics
!---
! (1) determine if today will rain or not
!---
!-- calculate monthly-average probability of rainy day
!---
    DO i=1,npoi
      pwet(i) = xinwet(i,imonth)/ndaypm(imonth)
    ENDDO
!---
    IF (.NOT.precip_exact) THEN
!-----
!---- (1.1) following Geng et al.
!-----
      IF (is_root_prc) THEN
         CALL random_number (rndnum)
      ENDIF
      CALL bcast(rndnum)
!-----
      DO i=1,npoi
        IF (xinprec(i,imonth) > 1.e-6) THEN
!---------
!-------- implement simple first-order Markov-chain precipitation
!-------- generator logic based on Geng et al. (1986),
!-------- Richardson and Wright (1984), and Richardson (1981)
!---------
!-------- basically, this allows for the probability of today being
!-------- a wet day (a day with measureable precipitation)
!-------- to be a function of what yesterday was (wet or dry)
!---------
!-------- the logic here is that it is more likely that a wet day
!-------- will follow another wet day -- allowing
!-------- for 'storm events' to persist
!---------
!-------- estimate the probability of a wet day after a dry day
!---------
          pwd = 0.75*pwet(i)
!---------
!-------- estimate the probability of a wet day after a wet day
!---------
          pww = 0.25+pwd
!---------
!-------- decide if today is a wet day or a dry day
!-------- using a random number
!---------
!-------- call random_number(rndnum) ! done before
!---------
          IF (iwet(i) == 0) then
            IF (rndnum <= pwd) iwet(i) = 1
          ELSE
            IF (rndnum > pww) iwet(i) = 0
          ENDIF
        ELSE
          iwet(i) = 0
        ENDIF
      ENDDO
    ELSE
!-----
!---- (1.2) preserving the monthly number of precip days
!----       and monthly precip
!-----
      DO i=1,npoi
        IF (ABS(xinwet(i,imonth)) < 32.) THEN
          IF (xinprec(i,imonth) > 1.e-6) THEN
            IF (   jour_precip(iday,imonth) &
 &              <= NINT(MAX(1.,xinwet(i,imonth))) ) THEN
              iwet(i) = 1
            ELSE
              iwet(i) = 0
            ENDIF
          ELSE
           iwet(i) = 0
          ENDIF
        ENDIF
      ENDDO
    ENDIF
!---
!-- (2) determine today's precipitation amount
!---
    IF (.not.precip_exact) THEN
       Warning_aa_ab(:)=.FALSE.
       Warning_iwet(:)=.FALSE.
!-----
!---- (2.1) following Geng et al.
!-----
      aa(:) = zero
      ab(:) = zero
      tr2(:)= zero
      tr1(:)= zero
      beta(:) = un
      DO i=1,npoi
!-------
!------ initialize daily precipitation to zero
!-------
        precip(i) = 0.0
!-------
        IF (xinprec(i,imonth) > 1.e-6) THEN
!---------
!-------- if it is going to rain today
!---------
          IF (iwet(i) == 1) THEN
!-----------
!---------- calculate average rainfall per wet day
!-----------
            rainpwd(i) = xinprec(i,imonth) &
 &                      *ndaypm(imonth)/MAX(0.1,xinwet(i,imonth))
!-----------
!---------- randomly select a daily rainfall amount
!---------- from a probability density function of rainfall
!----------
!---------- method i --
!-----------
!---------- use the following technique from Geng et al. and Richardson
!---------- to distribute rainfall probabilities
!-----------
!---------- pick a random rainfall amount
!---------- from a two-parameter gamma function distribution function
!-----------
!---------- estimate two parameters for gamma function
!---------- (following Geng et al.)
!-----------
            beta(i) = MAX(1.0,-2.16+1.83*rainpwd(i))
            alpha(i) = rainpwd(i)/beta(i)
!-----------
!---------- determine daily precipitation amount
!---------- from gamma distribution function
!---------- (following WGEN code of Richardson and Wright (1984))
!-----------
            IF (ABS(1.-alpha(i)) < 1.e-6) THEN
              alpha(i) = 1.e-6*(alpha(i)/ABS(alpha(i)))
            ENDIF
            aa(i) = 1.0/alpha(i)
            ab(i) = 1.0/(1.0-alpha(i))
!-----------
            IF ( (ABS(aa(i)) < 1.e-6) .OR. (ABS(ab(i)) < 1.e-6) ) THEN
               Warning_aa_ab(:)=.TRUE.
            ENDIF
            tr1(i) = exp(-18.42/aa(i))
            tr2(i) = exp(-18.42/ab(i))
          ENDIF
        ELSE
          IF (iwet(i) == 1) THEN
            Warning_iwet(i)=.TRUE.
          ENDIF
        ENDIF
      ENDDO

      DO i=1,npoi
         IF ( Warning_aa_ab(i) ) THEN
            WRITE(numout,*) ' ATTENTION, aa ou ab:'
            WRITE(numout,*) ' aa, ab = ',aa(i),ab(i)
            WRITE(numout,*) '   alpha, rainpwd, beta =', &
 &                       alpha(i),rainpwd(i),beta(i)
         ENDIF
         IF ( Warning_iwet(i) ) THEN
            WRITE(numout,*) ' ATTENTION, iwet = 1 alors que xinprec = 0)'
            WRITE(numout,*) '   xinprec, iwet = ',xinprec(i,imonth),iwet(i)
          ENDIF
      ENDDO
!-----
      where ( iwet(:) == 1 )
        not_ok(:) = 1
      elsewhere
        not_ok(:) = 0
      endwhere
!-----
      count_not_ok_g=0
      count_not_ok=SUM(not_ok(:))
      CALL reduce_sum(count_not_ok,count_not_ok_g)
      CALL bcast(count_not_ok_g)
!-
      z(:)=0.0
      DO WHILE (count_not_ok_g > 0)
        IF (is_root_prc) THEN
        CALL random_number (rn1)
        CALL random_number (rn2)
        ENDIF
        CALL bcast(rn1)
        CALL bcast(rn2)

        DO i=1,npoi
          IF ((iwet(i) == 1).AND.(not_ok(i) == 1)) then
            IF ( (rn1-tr1(i)) <= 0. ) THEN
              s1 = 0.0
            ELSE
              s1 = rn1**aa(i)
            ENDIF
!-----------
            IF ((rn2-tr2(i)) <= 0.) THEN
              s2 = 0.0
            ELSE
              s2 = rn2**ab(i)
            ENDIF
!-----------
            s12 = s1+s2
            IF ((s12-1.0) <= 0.) THEN
              z(i) = s1/s12
              not_ok(i) = 0
            ENDIF
          ENDIF
        ENDDO

        count_not_ok_g=0
        count_not_ok=SUM(not_ok(:))
        CALL reduce_sum(count_not_ok,count_not_ok_g)
        CALL bcast(count_not_ok_g)
      ENDDO
!-----
      IF (is_root_prc) THEN
         CALL random_number (rn3)
      ENDIF
      CALL bcast(rn3)
!      WRITE(*,*) mpi_rank,"rn3=",rn3
!-----
      DO i=1,npoi
        IF (iwet(i) == 1) then
          precip(i) = -z(i)*log(rn3)*beta(i)
        ENDIF
      ENDDO
!-----
!---- method ii --
!-----
!---- here we use a one-parameter Weibull distribution function
!---- following the analysis of Selker and Haith (1990)
!-----
!---- Selker, J.S. and D.A. Haith, 1990: Development and testing
!---- of single- parameter precipitation distributions,
!---- Water Resources Research, 11, 2733-2740.
!-----
!---- this technique seems to have a significant advantage over other
!---- means of generating rainfall distribution functions
!-----
!---- by calibrating the Weibull function to U.S. precipitation records,
!---- Selker and Haith were able to establish the following relationship
!-----
!---- the cumulative probability of rainfall intensity x is given as:
!-----
!---- f(x) = 1.0-exp(-(1.191 x / rainpwd)**0.75)
!-----
!---- where x       : rainfall intensity
!----       rainpwd : rainfall per wet day
!-----
!---- using transformation method, take uniform deviate and convert
!---- it to a random number weighted by the following Weibull function
!-----
!----    call random_number(rndnum)
!-----
!----    precip(i) = rainpwd / 1.191*(-log(1.0-rndnum))**1.333333
!-----
!---- bound daily precipitation to "REAListic" range
!-----
      DO i=1,npoi
        IF (iwet(i) == 1) THEN
!---------
!-------- lower end is determined by definition of a 'wet day'
!-------- (at least 0.25 mm of total precipitation)
!---------
!-------- upper end is to prevent ibis from blowing up
!---------
          precip(i) = MAX(precip(i),  0.25) ! min =   0.25 mm/day
          precip(i) = MIN(precip(i),150.00) ! max = 150.00 mm/day
        ENDIF
      ENDDO
    ELSE
!-----
!---- (2.2) preserving the monthly number of precip days
!----       and monthly precip
!-----
      DO i=1,npoi
!-------
!------ Correction Nathalie. C'est abs(xinwet) < 32 qu'il faut tester
!------ et non pas abs(xinprec(i,imonth)) < 32.
!-------
!------ IF ( (iwet(i) == 1).and.(abs(xinprec(i,imonth)) < 32.) ) THEN
        IF ( (iwet(i) == 1).and.(abs(xinwet(i,imonth)) < 32.) ) THEN
          precip(i) = xinprec(i,imonth)*REAL(ndaypm(imonth)) &
 &                   /REAL(NINT(MAX(1.,xinwet(i,imonth))))
        ELSE
          precip(i) = 0.
        ENDIF
      ENDDO
    ENDIF
!---
!-- (3) estimate expected minimum and maximum temperatures
!---
    DO i=1,npoi
!-----
!---- first determine the expected maximum and minimum temperatures
!---- (from climatological means) for this day of the year
!-----
!---- mean daily mean temperature (K)
      tdm = xint(i,it1w)+dt*(xint(i,it2w)-xint(i,it1w))+zero_t
!---- mean daily temperature range (K)
      trngm = xintrng(i,it1w)+dt*(xintrng(i,it2w)-xintrng(i,it1w))
!---- mean minimum and maximum temperatures
      tmaxm = tdm+0.56*trngm
      tminm = tdm-0.44*trngm
!-----
!---- modify maximum temperatures for wet and dry days
!-----
      IF (pwet(i) /= 0.0) THEN
        tmaxd = tmaxm+pwet(i)*omtmax*trngm
        tmaxw = tmaxd-        omtmax*trngm
      ELSE
        tmaxd = tmaxm
        tmaxw = tmaxm
      ENDIF
!-----
!---- set the 'expected' maximum and minimum temperatures for today
!-----
!---- note that the expected minimum temperatures are the same for
!---- both wet and dry days
!-----
      if (iwet(i) == 0)  tmaxe(i) = tmaxd
      if (iwet(i) == 1)  tmaxe(i) = tmaxw
!-----
      tmine(i) = tminm
!-----
!---- estimate variability in minimum and maximum temperatures
!-----
!---- tmaxs : standard deviation in maximum temperature (K)
!---- tmins : standard deviation in minimum temperature (K)
!-----
!---- Regression is based on analysis of 2-m air temperature data
!---- from the NCEP/NCAR reanalysis (1958-1997) for 294 land points
!---- over central North America
!---- (24N-52N, 130W-60W, 0.5-degree resolution):
!---- Daily max and min temperatures were calculated for each
!---- land point from daily mean temperature and temperature range.
!---- Anomalies were calculated
!---- by subtracting similar max and min temperatures calculated from
!---- monthly mean temperature and range (interpolated to daily).
!---- The 40 years of anomalies were then binned by month
!---- and the standard deviation calculated for each month.
!---- The 294 x 12 standard deviations were then regressed
!---- against the 3528 long-term monthly mean temperatures.
!-----
!---- note: the values are bound to be greater than 1.0 K
!----       (at the very least they must be bound
!----        so they don't go below zero)
!-----
      tmaxs(i) = MAX(1.0,-0.0713*(tdm-zero_t)+4.89)
      tmins(i) = MAX(1.0,-0.1280*(tdm-zero_t)+5.73)
    ENDDO
!---
!-- (4) estimate expected cloud cover
!---
!---
!-- the formulation of dry and wet cloud cover has been
!-- derived from the weather generator used in the epic crop model
!---
    DO i=1,npoi
!-----
!---- cloudm : mean cloud cover for today
!---- cloudd : dry day cloud cover
!---- cloudw : dry day cloud cover
!---- cloude : expected cloud cover today
!-----
!---- mean cloud cover (%)
!-----
      cloudm = xincld(i,it1w)+dt*(xincld(i,it2w)-xincld(i,it1w))
!-----
!---- convert from percent to fraction
!-----
      cloudm = cloudm/100.0
!-----
!---- adjust cloud cover depending on dry day / rainy day
!---- following logic of the EPIC weather generator code
!-----
      IF (pwet(i) /= 0.0) THEN
        cloudd = (cloudm-pwet(i)*omcloud)/(1.0-pwet(i)*omcloud)
        cloudd = MIN(1.0,MAX(0.0,cloudd))
        cloudw = (cloudm-(1.0-pwet(i))*cloudd)/pwet(i)
      ELSE
        cloudd = cloudm
        cloudw = cloudm
      ENDIF
      IF (iwet(i) == 0)  cloude(i) = cloudd
      IF (iwet(i) == 1)  cloude(i) = cloudw
!-----
!---- estimate variability in cloud cover for wet and dry days
!---- following numbers proposed by Richardson
!-----
!---- clouds : standard deviation of cloud fraction
!-----
      IF (iwet(i) == 0)  clouds(i) = 0.24*cloude(i)
      IF (iwet(i) == 1)  clouds(i) = 0.48*cloude(i)
    ENDDO
!---
! (5) determine today's temperatures and cloud cover using
!     first-order serial autoregressive technique
!---
!-- use the Richardson (1981) weather generator approach to simulate the
!-- daily values of minimum / maximum temperature and cloud cover
!---
!-- following the implementation of the Richardson WGEN weather
!-- generator used in the EPIC crop model
!---
!-- this approach uses a multivariate generator, which assumes that the
!-- perturbation of minimum / maximum temperature and cloud cover are
!-- normally distributed and that the serial correlation of each
!-- variable may be described by a first-order autoregressive model
!---
!-- generate standard deviates for weather generator
!---
    DO j=1,3
      ee(j) = 9999.
      DO WHILE (ee(j) > 2.5)
        IF (is_root_prc) THEN
           CALL random_number (rn1)
           CALL random_number (rn2)
        ENDIF
        CALL bcast(rn1)
        CALL bcast(rn2)
        ee(j) = SQRT(-2.0*LOG(rn1))*COS(6.283185*rn2)
      ENDDO
    ENDDO
!---
!-- zero out vectors
!---
    r(1:3)  = 0.0
    rr(1:npoi,1:3) = 0.0
!---
!-- update working vectors
!---
    do j=1,3
      do k=1,3
        r(j) = r(j)+b(j,k)*ee(j)
      enddo
    enddo
!---
    do j=1,3
      do k=1,3
        do i=1,npoi
          rr(i,j) = rr(i,j)+a(j,k)*xstore(i,k)
        enddo
      enddo
    enddo
!---
!-- solve for x() perturbation vector and save current vector
!-- into the xim1() storage vector (saved for each point)
!---
    do j=1,3
      do i=1,npoi
        xstore(i,j) = r(j)+rr(i,j)
      enddo
    enddo
!---
!-- determine today's minimum and maximum temperature
!--
    do i=1,npoi
      tmax(i) = tmaxe(i)+tmaxs(i)*xstore(i,1)
      tmin(i) = tmine(i)+tmins(i)*xstore(i,2)
!-----
!---- if tmin > tmax, then switch the two around
!-----
      if (tmin(i) > tmax(i)) then
        tdum    = tmax(i)
        tmax(i) = tmin(i)
        tmin(i) = tdum
      ENDIF
!---- daily average temperature
      td(i) = 0.44*tmax(i)+0.56*tmin(i)
!---- determine today's cloud cover
      cloud(i) = cloude(i)+clouds(i)*xstore(i,3)
!---- constrain cloud cover to be between 0 and 100%
      cloud(i) = MAX(0.0,MIN(1.0,cloud(i)))
    enddo
!---
!-- (6) estimate today's surface atmospheric pressure
!---
    do i=1,npoi
!-----
!---- simply a function of the daily average temperature
!---- and topographic height -- nothing fancy here
!-----
      psurf(i) = 101325.0*(td(i)/(td(i)+0.0065*xintopo(i)))**rwork
    enddo
!---
!-- (7) estimate today's relative humidity
!---
    IF (is_root_prc) THEN
       CALL random_number (rn)
    ENDIF
    CALL bcast(rn)
!---
    CALL weathgen_qsat (npoi,td,psurf,qsattd)
!---
    do i=1,npoi
!-----
!---- the formulation of dry and wet relative humidities has been
!---- derived from the weather generator used in the epic crop model
!-----
!---- qdm : mean relative humidity
!---- qdd : dry day relative humidity
!---- qdw : rainy day relative humidity
!---- qde : expected relative humidity (based on wet/dry decision)
!-----
!---- mean relative humidity (%)
      qdm(i) = xinq(i,it1w)+dt*(xinq(i,it2w)-xinq(i,it1w))
!---- convert from percent to fraction
      qdm(i) = qdm(i)/100.0
!-----
!---- adjust humidity depending on dry day / rainy day
!---- following logic of the EPIC weather generator code
!-----
      if (pwet(i) /= 0.0) then
        qdd(i) = (qdm(i)-pwet(i)*omqd)/(1.0-pwet(i)*omqd)
        if (qdd(i) < 0.2) then
          qdd(i) = 0.2
          if (qdd(i) > qdm(i)) qdm(i) = qdd(i)
        ENDIF
        qdd(i) = MIN(1.0,qdd(i))
        qdw(i) = (qdm(i)-(1.0-pwet(i))*qdd(i))/pwet(i)
      ELSE
        qdd(i) = qdm(i)
        qdw(i) = qdm(i)
      ENDIF
!-----
      if (iwet(i) == 0)  qde(i) = qdd(i)
      if (iwet(i) == 1)  qde(i) = qdw(i)
!-----
!---- estimate lower and upper bounds of humidity distribution function
!---- following logic of the EPIC weather generator code
!-----
      xx = exp(qde(i))
      qdup(i)  = qde(i)+(1.0-qde(i))*xx/e
      qdlow(i) = qde(i)*(1.0-1./xx)
!-----
!---- randomlly select humidity from triangular distribution function
!---- following logic of the EPIC weather generator code
!-----
!---- call random_number(rn) ! GK done before
!-----
      y  = 2.0/(qdup(i)-qdlow(i))
!-----
      b3 = qde(i)-qdlow(i)
      b2 = qdup(i)-qde(i)
      b1 = rn/y
!-----
      x1 = y*b3/2.0
!-----
      if (rn > x1) then
        qd(i) = qdup(i)-sqrt (b2*b2-2.0*b2*(b1-0.5*b3))
      ELSE
        qd(i) = qdlow(i)+sqrt (2.0*b1*b3)
      ENDIF
!-----
!---- adjust daily humidity to conserve monthly mean values
!-----
!---- note that this adjustment sometimes gives rise to humidity
!---- values greater than 1.0 -- which is corrected below
!-----
      amn = (qdup(i)+qde(i)+qdlow(i))/3.0
      qd(i) = qd(i)*qde(i)/amn
!-----
!---- constrain daily average relative humidity
!-----
      qd(i) = MAX(0.30,MIN(qd(i),0.99))
!-----
!---- convert from relative humidity to specific humidity at
!---- daily mean temperature
!-----
      qd(i) = qd(i)*qsattd(i)
    enddo
!---
!-- (8) estimate today's daily average wind speed
!---
    IF (is_root_prc) THEN
        CALL random_number (rn4)
    ENDIF
    CALL bcast(rn4)

    do i=1,npoi
!-----
!---- first estimate the expected daily average wind speed (from monthly
!---- means)
!-----
      eud = xinwind(i,it1w)+dt*(xinwind(i,it2w)-xinwind(i,it1w))
!-----
!---- following logic of the EPIC weather generator
!---- select random wind speed following this equation
!-----
!---- call random_number(rn4)
!-----
      ud(i) = 1.13989*eud*(-log(rn4))**0.30
!---- constrain daily wind speeds to be between 2.5 and 10.0 m/sec
      ud(i) = MAX(2.5,MIN(ud(i),10.0))
    ENDDO
  ELSE
!---
!-- use REAL daily climate data
!---
    DO i=1,npoi
!-----
!---- use basic daily climate data, converting units
!-----
!---- daily total precipitation
      precip(i) = xinprec(i,imonth)
!---- daily average temperatures
      td(i) = xint(i,imonth)+zero_t
      trngm = MIN(44.0,xintrng(i,imonth))
!-----
      tmax(i) = td(i)+0.56*trngm
      tmin(i) = td(i)-0.44*trngm
!---- daily average cloud cover
      cloud(i) = xincld(i,imonth)*0.01
!---- daily average specific humidity
      qd(i) = xinq(i,imonth)
!---- daily average wind speed
      ud(i) = xinwind(i,imonth)
!-----
!---- compute surface atmospheric pressure
!-----
      psurf(i) = 101325.0*(td(i)/(td(i)+0.0065*xintopo(i)))**rwork
    ENDDO
  ENDIF
!-------------------
END SUBROUTINE daily
!-
!===
!-

!! ================================================================================================================================
!! SUBROUTINE   : diurnal
!!
!>\BRIEF          This routine interpolates values from the!* 
!! subroutine daily into instantaneous values for simulating a diurnal cycle.
!!
!! DESCRIPTION  :  
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S):  ::fira, ::solad, ::ua, ::ta, ::qa, ::raina, ::snowa, ::rh.
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE diurnal &
& (npoi, nband, time, jday, plens, startp, endp, latitude, &
&  cloud, tmax, tmin, precip, qd, ud, psurf, &
&  fira, solad, solai, ua, ta, qa, raina, snowa, rh)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-

  !! 0. Parameter and variables declaration

  REAL,PARAMETER :: stef = 5.67051E-8  !! Stefan-Boltzman constant @tex $(W.m^{-2}.K^{-4})$ @endtex 
  REAL,PARAMETER :: pi = 3.1415927     !! pi
  REAL,PARAMETER :: pir = pi/180.      !! conversion factor degrees to radians

  !! 0.1 Input variables

  INTEGER,INTENT(IN) :: npoi    !! number of grid points (unitless)
  INTEGER,INTENT(IN) :: nband   !! number of visible bands (unitless)
  REAL,INTENT(IN) :: time
  INTEGER, INTENT(IN) :: jday
  REAL,INTENT(IN) :: plens,startp,endp
  REAL,INTENT(IN) :: latitude(npoi)   !! latitude in degrees
  REAL,INTENT(IN) :: cloud(npoi)      !! cloud fraction (0-1, unitless)
  REAL,INTENT(IN) :: tmax(npoi)       !! maximum daily temperature (K)
  REAL,INTENT(IN) :: tmin(npoi)       !! maximum daily temperature (K)
  REAL,INTENT(IN) :: precip(npoi)     !! daily precitation @tex $(mm.day^{-1})$ @endtex 
  REAL,INTENT(IN) :: qd(npoi)         !! daily average specific humidity (kg_h2o/kg_air)
  REAL,INTENT(IN) :: ud(npoi)         !! daily average wind speed @tex $(m.s^{-1})$ @endtex 
  REAL,INTENT(IN) :: psurf(npoi)      !! surface pressure (Pa)

  !! 0.2 Output variables

  REAL,INTENT(OUT) :: fira(npoi)          !! incoming ir flux  @tex $(W.m^{-2})$ @endtex
  REAL,INTENT(OUT) :: solad(npoi,nband)   !! direct downward solar flux @tex $(W.m^{-2})$ @endtex
  REAL,INTENT(OUT) :: solai(npoi,nband)   !! diffuse downward solar flux @tex $(W.m^{-2})$ @endtex
  REAL,INTENT(OUT) :: ua(npoi)            !! wind speed @tex $(m.s^{-1})$ @endtex 
  REAL,INTENT(OUT) :: ta(npoi)            !! air temperature (K)
  REAL,INTENT(OUT) :: qa(npoi)            !! specific humidity (kg_h2o/kg_air)
  REAL,INTENT(OUT) :: raina(npoi)         !! rainfall rate @tex $(mm.day^{-1})$ @endtex 
  REAL,INTENT(OUT) :: snowa(npoi)         !! snowfall rate @tex $(mm.day^{-1})$ @endtex 
  REAL,INTENT(OUT) :: rh(npoi)            !! relative humidity (0-1, unitless)

  !! 0.4 Local variables

  REAL,SAVE      :: step
  REAL :: xl,so,xllp,xee,xse
  REAL :: xlam,dlamm,anm,ranm,ranv,anv,tls,rlam
  REAL :: sd,cd,deltar,delta,Dis_ST,ddt
!-
  REAL :: coszen(npoi)                      !! cosine of solar zenith angle
  REAL :: rtime
  REAL :: orbit,angle,xdecl,xlat
  REAL :: sw,frac,gamma,qmin,qmax,qsa,emb,ea,ec,dtair,dtcloud,rn
  REAL :: trans(npoi), fdiffuse(npoi), qsatta(npoi), qsattmin(npoi)
  INTEGER :: i,ib
  INTEGER,SAVE :: npoi0
  LOGICAL,SAVE :: firstcall = .TRUE.

!---------------------------------------------------------------------
! GK240100
  IF (firstcall) THEN
     IF ( TRIM(calendar_str) .EQ. 'gregorian' ) THEN  
        step = 1.0
     ELSE
        step = one_year/365.2425
     ENDIF
    firstcall = .FALSE.
    npoi0 = npoi
  ELSE IF (npoi /= npoi0) THEN
    WRITE(numout,*) 'Domain size old, new: ',npoi0,npoi
    STOP 'WG Diurnal: Problem: Domain has changed since last call'
  ENDIF
!-
! calendar and orbital calculations
!-
! calculate time in hours
  rtime = time/3600.0
!-
! calculate the earth's orbital angle (around the sun) in radians
  orbit = 2.0*pi*REAL(jday)/365.2425
!-
! calculate the solar hour angle in radians
  angle = 2.0*pi*(rtime-12.0)/24.0
!-
! calculate the current solar declination angle
! ref: global physical climatology, hartmann, appendix a
!
! xdecl = 0.006918 &
!        -0.399912*cos(orbit) &
!        +0.070257*sin(orbit) &
!        -0.006758*cos(2.0*orbit) &
!        +0.000907*sin(2.0*orbit) &
!        -0.002697*cos(3.0*orbit) &
!        +0.001480*sin(3.0*orbit)
!
! calculate the effective solar constant,
! including effects of eccentricity
! ref: global physical climatology, hartmann, appendix a
!
! sw = 1370.*( 1.000110 &
!             +0.034221*cos(orbit) &
!             +0.001280*sin(orbit) &
!             +0.000719*cos(2.0*orbit) &
!             +0.000077*sin(2.0*orbit))
!
! correction Marie-France Loutre
!
!    orbital parameters
!
!    ecc = 0.016724
!    perh = 102.04
!    xob = 23.446
!-
  xl = perh+180.0
! so : sinus of obliquity
  so = sin(xob*pir)
!-
  xllp = xl*pir
  xee  = ecc*ecc
  xse  = sqrt(1.0d0-xee)
! xlam : true long. sun for mean long. = 0
  xlam = (ecc/2.0+ecc*xee/8.0d0)*(1.0+xse)*sin(xllp)-xee/4.0 &
 &      *(0.5+xse)*sin(2.0*xllp)+ecc*xee/8.0*(1.0/3.0+xse) &
 &      *sin(3.0*xllp)
  xlam  =2.0d0*xlam/pir
! dlamm : mean long. sun for ma-ja
  dlamm =xlam+(jday-79)*step
  anm   = dlamm-xl
  ranm  = anm*pir
  xee   = xee*ecc
! ranv : anomalie vraie    (radian)
  ranv = ranm+(2.0*ecc-xee/4.0)*sin(ranm)+5.0/4.0*ecc*ecc &
 &      *sin(2.0*ranm)+13.0/12.0*xee*sin(3.0*ranm)
! anv  : anomalie vraie    (degrees)
  anv  = ranv/pir
! tls  : longitude vraie   (degrees)
  tls  = anv+xl
! rlam : longitude vraie   (radian)
  rlam = tls*pir
! sd and cd: cosinus and sinus of solar declination angle (delta)
! sinus delta = sin (obl)*sin(lambda) with lambda = real longitude
! (Phd. thesis of Marie-France Loutre, ASTR-UCL, Belgium, 1993)
  sd    = so*sin(rlam)
  cd    = sqrt(1.0d0-sd*sd)
! delta : Solar Declination (degrees, angle sun at equator)
  deltar = atan(sd/cd)
  delta  = deltar/pir
!-
! Eccentricity Effect
!-
  Dis_ST  =(1.0-ecc*ecc)/(1.0+ecc*cos(ranv))
! ddt :    1 / normalized earth's sun distance
  ddt = 1.0/Dis_ST
!-
! Insolation normal to the atmosphere (W/m2)
!-
  sw  = ddt *ddt *1370.d0
!-
  xdecl = deltar
!-
! solar calculations
!-
  do i=1,npoi
!---
!-- calculate the latitude in radians
!---
    xlat = latitude(i)*pi/180.0
!---
!-- calculate the cosine of the solar zenith angle
!---
    coszen(i) = MAX(0.0, (sin(xlat)*sin(xdecl) &
 &                      + cos(xlat)*cos(xdecl)*cos(angle)))
!---
!-- calculate the solar transmission through the atmosphere
!-- using simple linear function of tranmission and cloud cover
!---
!-- note that the 'cloud cover' data is typically obtained from
!-- sunshine hours -- not direct cloud observations
!---
!-- where, cloud cover = 1 - sunshine fraction
!---
!-- different authors present different values for the slope and
!-- intercept terms of this equation
!---
!-- Friend, A: Parameterization of a global daily weather generator
!-- for terrestrial ecosystem and biogeochemical modelling,
!-- Ecological Modelling
!---
!-- Spitters et al., 1986: Separating the diffuse and direct component
!-- of global radiation and its implications for modeling canopy
!-- photosynthesis, Part I: Components of incoming radiation,
!-- Agricultural and Forest Meteorology, 38, 217-229.
!---
!-- A. Friend       : trans = 0.251+0.509*(1.0-cloud(i))
!-- Spitters et al. : trans = 0.200+0.560*(1.0-cloud(i))
!---
!-- we are using the values from A. Friend
!---
    trans(i) = 0.251+0.509*(1.0-cloud(i))
!---
!-- calculate the fraction of indirect (diffuse) solar radiation
!-- based upon the cloud cover
!---
!-- note that these relationships typically are measured for either
!-- monthly or daily timescales, and may not be exactly appropriate
!-- for hourly calculations -- however, in ibis, cloud cover is fixed
!-- through the entire day so it may not make much difference
!---
!-- method i --
!---
!-- we use a simple empirical relationships
!-- from Nikolov and Zeller (1992)
!---
!-- Nikolov, N. and K.F. Zeller, 1992:  A solar radiation algorithm
!-- for ecosystem dynamics models, Ecological Modelling, 61, 149-168.
!---
    fdiffuse(i) = 1.0045+trans(i)*( 0.0435+trans(i) &
 &                                 *(-3.5227+trans(i)*2.6313))
!---
    IF (trans(i) > 0.75) fdiffuse(i) = 0.166
!---
!-- method ii --
!---
!-- another method was suggested by Spitters et al. (1986) based on
!-- long-term data from the Netherlands
!--
!-- Spitters et al., 1986: Separating the diffuse and direct component
!-- of global radiation and its implications for modeling canopy
!-- photosynthesis, Part I: Components of incoming radiation,
!-- Agricultural and Forest Meteorology, 38, 217-229.
!---
!--       if ((trans == 0.00).and.(trans < 0.07)) then
!--         fdiffuse = 1.0
!--       else if ((trans >= 0.07).and.(trans < 0.35)) then
!--         fdiffuse = 1.0-2.3*(trans-0.07)**2
!--       else if ((trans >= 0.35).and.(trans < 0.75)) then
!--         fdiffuse = 1.33-1.46*trans
!--       ELSE
!--         fdiffuse = 0.23
!--       endif
!---
  ENDDO
!-
! do for each waveband
!-
  DO ib=1,nband
!---
!-- calculate the fraction in each waveband
!---
!-- GK010200
    IF (nband == 2) then
      frac = 0.46+0.08*REAL(ib-1)
    ELSE
      frac = 1./REAL(nband)
    ENDIF
!---
    DO i=1,npoi
!-----
!---- calculate the direct and indirect solar radiation
!-----
      solad(i,ib) = sw*coszen(i)*frac*trans(i)*(1.-fdiffuse(i))
      solai(i,ib) = sw*coszen(i)*frac*trans(i)*fdiffuse(i)
    ENDDO
  ENDDO
!-
! temperature calculations
!-
!---
!-- assign hourly temperatures using tmax and tmin
!-- following Environmental Biophysics, by Campbell and Norman, p.23
!---
!-- this function fits a fourier series to the diurnal temperature cycle
!-- note that the maximum temperature occurs at 2:00 pm local solar time
!---
!-- note that the daily mean value of gamma is 0.44,
!-- so td = 0.44*tmax+0.56*tmin,  instead of
!--    td = 0.50*tmax+0.50*tmin
!---
  gamma = 0.44-0.46*SIN(    pi/12.0*rtime+0.9) &
              +0.11*SIN(2.0*pi/12.0*rtime+0.9)
  DO i=1,npoi
    ta(i) = tmax(i)*gamma+tmin(i)*(1.0-gamma)
  ENDDO
!-
! humidity calculations
!-
  CALL weathgen_qsat (npoi,tmin,psurf,qsattmin)
  CALL weathgen_qsat (npoi,ta,psurf,qsatta)
!-
  DO i=1,npoi
!---
!-- adjust specific humidity against daily minimum temperatures
!---
!-- To do this, qa is written as an approximate sine function
!-- (same as ta) to preserve the daily mean specific humidity,
!-- while also preventing rh from exceeding 99% at night
!---
!-- Note that the daily mean RH is *not* preserved, and therefore the
!-- output RH will be slightly different from what was read in.
!---
!-- first adjust so that maximum RH cannot exceed 99% at night
!---
    qmin = MIN(qd(i),0.99*qsattmin(i))
    qmax = (qd(i)-0.56*qmin)/0.44
!---
!-- if needed, adjust again to 99% at other times of the day (in which
!-- case the daily mean *specific* humidity is also not preserved)
!---
    qsa  = 0.99*qsatta(i)
!---
!-- calculate the hourly specific humidity, using the above adjustments
!---
    qa(i) = MIN(qsa,qmax*gamma+qmin*(1.0-gamma))
!---
!-- calculate the hourly relative humidity
!--
    rh(i) = 100.*qa(i)/qsatta(i)
  ENDDO
!-
! wind speed calculations
!-
  IF (is_root_prc) THEN
     CALL random_number (rn)
  ENDIF
  CALL bcast(rn)
!-
  DO i=1,npoi
!---
!-- following logic of the EPIC weather generator
!-- select random wind speed following this equation
!---
!-- call random_number(rn) ! done before
!---
    ua(i) = 1.13989*ud(i)*(-log(rn))**0.30
!---
!-- fix wind speeds to always be above 2.5 m/sec and below 10.0 m/sec
!---
    ua(i) = MAX(2.5,MIN(10.0,ua(i)))
  ENDDO
!-
! ir flux calculations
!-
  DO i=1,npoi
!---
!-- clear-sky emissivity as a function of water vapor pressure
!-- and atmospheric temperature
!---
!-- calculate the ir emissivity of the clear sky
!-- using equation from idso (1981) water resources res., 17, 295-304
!---
    emb = 0.01*(psurf(i)*qa(i)/(0.622+qa(i)))
    ea  = 0.70+5.95e-5*emb*EXP(1500.0/ta(i))
!---
!-- assign the ir emissivity of clouds
!-- (assume they are ~black in the ir)
!---
    ec = 0.950
!---
!-- assign the temperature difference of emissions (air+cloud) from
!-- the surface air temperature
!---
    dtair   = 0.0
    dtcloud = 0.0
!---
!-- total downward ir is equal to the sum of:
!---
!-- (1) clear sky contribution to downward ir radiation flux
!-- (2) cloud contribution to downward ir radiation flux
!---
    fira(i) = (1.-cloud(i))*ea*stef*(ta(i)-dtair)**4 &
 &           +cloud(i)*ec*stef*(ta(i)-dtcloud)**4
  ENDDO
!-
! snow and rain calculations
!-
  DO i=1,npoi
!---
!-- reset snow and rain to zero
!---
    snowa(i) = 0.0
    raina(i) = 0.0
!---
!-- if precipitation event then calculate
!---
    IF (time >= startp .and. time < endp) then
!-----
!---- for rain / snow partitioning, make it all rain if
!---- ta > 2.5 C and all snow if ta <= 2.5 C
!-----
!---- reference:
!-----
!---- Auer, A. H., 1974: The rain versus snow threshold temperatures,
!---- Weatherwise, 27, 67.
!-----
      IF (ta(i)-273.15 > 2.5) then
        raina(i) = precip(i)/plens
      ELSE
        snowa(i) = precip(i)/plens
      ENDIF
    ENDIF
  ENDDO
!---------------------
END SUBROUTINE diurnal
!-
!===
!-

!! ================================================================================================================================
!! SUBROUTINE   : weathgen_main
!!
!>\BRIEF        This subroutine manages the calls to the different subroutines of the weather module. 
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::tair, ::pb, ::qair, ::swdown, ::rainf, ::snowf, ::u, ::v, ::lwdown
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE weathgen_main &
& (itau, istp, filename, force_id, iim, jjm, nband, &
&  rest_id, lrstread, lrstwrite, &
&  limit_west, limit_east, limit_north, limit_south, &
&  zonal_res, merid_res, lon, lat, date0, dt_force, &
&  kindex, nbindex, &
&  swdown, rainf, snowf, tair, u, v, qair, pb, lwdown)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-

  !! 0. Variables and parameters declaration

  !! 0.1 Input variables

  INTEGER,INTENT(IN)                  :: itau,istp
  CHARACTER(LEN=*),INTENT(IN)         :: filename
  INTEGER,INTENT(IN)                  :: force_id
  INTEGER,INTENT(IN)                  :: iim,jjm
  INTEGER,INTENT(IN)                  :: nband
  INTEGER,INTENT(IN)                  :: rest_id
  LOGICAL,INTENT(IN)                  :: lrstread, lrstwrite
  REAL,INTENT(IN)                     :: limit_west,limit_east
  REAL,INTENT(IN)                     :: limit_north,limit_south
  REAL,INTENT(IN)                     :: zonal_res,merid_res
  REAL,INTENT(IN)                     :: date0,dt_force
  REAL,DIMENSION(iim,jjm),INTENT(IN)  :: lon,lat  !! longitude, latitude

  !! 0.2 Output variables 

  REAL,DIMENSION(iim,jjm),INTENT(OUT) :: &
 &  tair,pb,qair,swdown,rainf,snowf, u,v,lwdown

  !! 0.3 Mofified variables

  INTEGER,INTENT(INOUT)                    :: nbindex
  INTEGER,DIMENSION(iim*jjm),INTENT(INOUT) :: kindex

  !! 0.4 Local variables

  REAL, ALLOCATABLE, DIMENSION(:,:) :: &
 &  tair_g,pb_g,qair_g,swdown_g,rainf_g,snowf_g, u_g,v_g,lwdown_g
  REAL,DIMENSION(nbindex) :: &
 &  tairl,pbl,qairl,swdownl,rainfl,snowfl, ul,vl,lwdownl
  INTEGER :: i,j,ij
!---------------------------------------------------------------------
  IF (lrstread) THEN
    CALL weathgen_begin ( &
     & dt_force,itau, date0, &
     & rest_id,iim,jjm, &
     & lon,lat,tair,pb,qair,kindex,nbindex)
        RETURN
  ELSE
    CALL weathgen_get &
 &   (itau, date0, dt_force, nbindex, nband, lat_land, &
 &    swdownl, rainfl, snowfl, tairl, ul, vl, qairl, pbl, lwdownl)

    tair(:,:)=val_exp
    qair(:,:)=val_exp
    pb(:,:)=val_exp
!!$    tair(:,:)=280.
!!$    qair(:,:)=1.E-03
!!$    pb(:,:)=101325
    DO ij=1,nbindex
       j = (((kindex(ij)-1)/iim) + 1)
       i = (kindex(ij) - (j-1)*iim)
       !
       swdown(i,j)=swdownl(ij)
       rainf(i,j)=rainfl(ij)
       snowf(i,j)=snowfl(ij)
       tair(i,j)=tairl(ij)
       u(i,j)=ul(ij)
       v(i,j)=vl(ij)
       qair(i,j)=qairl(ij)
       pb(i,j)=pbl(ij)
       lwdown(i,j)=lwdownl(ij)
    ENDDO
!---
    IF (dump_weather) THEN
      ALLOCATE(tair_g(iim_g,jjm_g))
      ALLOCATE(pb_g(iim_g,jjm_g))
      ALLOCATE(qair_g(iim_g,jjm_g))
      ALLOCATE(swdown_g(iim_g,jjm_g))
      ALLOCATE(rainf_g(iim_g,jjm_g))
      ALLOCATE(snowf_g(iim_g,jjm_g))
      ALLOCATE(u_g(iim_g,jjm_g))
      ALLOCATE(v_g(iim_g,jjm_g))
      ALLOCATE(lwdown_g(iim_g,jjm_g))

      CALL gather2D(tair, tair_g)
      CALL gather2D(pb, pb_g)
      CALL gather2D(qair, qair_g)
      CALL gather2D(swdown, swdown_g)
      CALL gather2D(rainf, rainf_g)
      CALL gather2D(snowf, snowf_g)
      CALL gather2D(u, u_g)
      CALL gather2D(v, v_g)
      CALL gather2D(lwdown, lwdown_g)
      IF (is_root_prc) THEN
         CALL weathgen_dump &
 &           (itau, dt_force, iim_g, jjm_g, nbp_glo, index_g, lrstwrite, &
 &            swdown_g, rainf_g, snowf_g, tair_g, u_g, v_g, qair_g, pb_g, lwdown_g)
      ENDIF
    ENDIF
!---
    IF (lrstwrite) THEN
      CALL weathgen_restwrite (rest_id,istp,iim,jjm,nbindex,kindex)
    ENDIF
  ENDIF
!---------------------------
END SUBROUTINE weathgen_main
!-
!===
!-

!! ================================================================================================================================
!! SUBROUTINE   : weathgen_init 
!!
!>\BRIEF         This subroutine initializes weather variables.
!!
!! DESCRIPTION  : None 
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S):  ::kindex, ::nbindex
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE weathgen_init &
     & (filename,dt_force,force_id,iim,jjm, &
     &  zonal_res,merid_res,lon,lat,kindex,nbindex)
!!$,iind,jind)
  !---------------------------------------------------------------------
  IMPLICIT NONE
  !-

  !! 0. Variables and parameters declaration

  REAL,PARAMETER  :: fcrit = .5       !! Minimum land fraction on original grid

  !! 0.1 Input variables

  CHARACTER(LEN=*),INTENT(IN)         :: filename
  REAL,INTENT(IN)                     :: dt_force
  INTEGER,INTENT(IN)                  :: iim, jjm
  REAL,INTENT(IN)                     :: zonal_res,merid_res
  REAL,DIMENSION(iim,jjm),INTENT(IN)  :: lon,lat

  !! 0.2 Output variables  

  INTEGER,INTENT(OUT)                 :: nbindex
  INTEGER,DIMENSION(iim*jjm),INTENT(OUT)  :: kindex
!!$  INTEGER,DIMENSION(iim),INTENT(OUT) :: iind
!!$  INTEGER,DIMENSION(jjm),INTENT(OUT) :: jind

  !! 0.3 Modified variables 

  INTEGER,INTENT(INOUT)               :: force_id

  !! 0.4 Local variables

  REAL,DIMENSION(:),ALLOCATABLE         :: lon_file, lon_temp
  REAL,DIMENSION(:,:),ALLOCATABLE       :: nav_lon, nav_lat
  REAL,DIMENSION(:),ALLOCATABLE         :: lat_file, lat_temp
  REAL,DIMENSION(:,:),ALLOCATABLE       :: lsm_file
  REAL     :: xextent_file, yextent_file, xres_file, yres_file
  INTEGER  :: i, j, plotstep
  REAL,DIMENSION(iim,jjm)               :: mask
  CHARACTER(LEN=1),DIMENSION(0:1)       :: maskchar
  CHARACTER(LEN=30)                     :: var_name
  REAL     :: x_cut
  REAL,DIMENSION(iim) :: tmp_lon
  REAL,DIMENSION(jjm) :: tmp_lat

!_ ================================================================================================================================

  !-
  ! 0. messages, initialisations
  !-
  WRITE(numout,*) &
       &  'weathgen_init: Minimum land fraction on original grid =',fcrit
  CALL ioget_calendar(calendar_str)
  !-
  ! 1. on lit les longitudes et latitudes de la grille de depart
  !    ainsi que le masque terre-ocean
  !-
  CALL flinclo(force_id)
  CALL flininfo (filename,iim_file,jjm_file,llm_file,ttm_file,force_id)
  !-
  ALLOC_ERR=-1
  ALLOCATE(nav_lon(iim_file,jjm_file), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
     WRITE(numout,*) "ERROR IN ALLOCATION of nav_lon : ",ALLOC_ERR
     STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(nav_lat(iim_file,jjm_file), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
     WRITE(numout,*) "ERROR IN ALLOCATION of nav_lat : ",ALLOC_ERR
     STOP 
  ENDIF
  !-
  var_name='nav_lon'
  CALL flinget (force_id,var_name,iim_file,jjm_file,0,0,1,1,nav_lon)
  var_name='nav_lat'
  CALL flinget (force_id,var_name,iim_file,jjm_file,0,0,1,1,nav_lat)
  !-

  ALLOC_ERR=-1
  ALLOCATE(lon_file(iim_file), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
     WRITE(numout,*) "ERROR IN ALLOCATION of lon_file : ",ALLOC_ERR
     STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(lat_file(jjm_file), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
     WRITE(numout,*) "ERROR IN ALLOCATION of lat_file : ",ALLOC_ERR
     STOP 
  ENDIF
  !-
  DO i=1,iim_file
     lon_file(i) = nav_lon(i,1)
  ENDDO
  DO j=1,jjm_file
     lat_file(j) = nav_lat(1,j)
  ENDDO
!-

  ALLOC_ERR=-1
  ALLOCATE(lsm_file(iim_file,jjm_file), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of lsm_file : ",ALLOC_ERR
    STOP 
  ENDIF
!-
  var_name='lsmera'
  CALL flinget (force_id,var_name,iim_file,jjm_file,0,0,1,1,lsm_file)
!-
! 2. La resolution du modele ne doit pas etre superieure
!    a celle de la grille de depart
!-
  xextent_file = lon_file(iim_file)-lon_file(1)
  yextent_file = lat_file(1)-lat_file(jjm_file)
  xres_file = xextent_file/REAL(iim_file-1)
  yres_file = yextent_file/REAL(jjm_file-1)
!-
  IF (xres_file > zonal_res) THEN
    WRITE(numout,*) 'Zonal resolution of model grid too fine.'
    WRITE(numout,*) 'Resolution of original data (deg): ', xres_file
    STOP
  ENDIF
!-
  IF (yres_file > merid_res) THEN
    WRITE(numout,*) 'Meridional resolution of model grid too fine.'
    WRITE(numout,*) 'Resolution of original data (deg): ', yres_file
    STOP
  ENDIF
!-
! 3. On verifie la coherence des coordonnees de depart et d'arrivee.
!    Sinon, il faut deplacer une partie du monde (rien que ca).
!-
  i_cut = 0
!-

  ALLOC_ERR=-1
  ALLOCATE(lon_temp(iim_file), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of lon_temp : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(lat_temp(jjm_file), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of lat_temp : ",ALLOC_ERR
    STOP 
  ENDIF
!-
  IF (lon(iim,1) > lon_file(iim_file)) THEN
!-- A l'Est de la region d'interet, il n'y a pas de donnees
!-- le bout a l'Ouest de la region d'interet est deplace de 360 deg
!-- vers l'Est
    x_cut = lon(1,1)
    DO i=1,iim_file
      IF (lon_file(i) < x_cut) i_cut = i
    ENDDO
    IF ((i_cut < 1).OR.(i_cut >= iim_file)) THEN
        STOP 'Cannot find longitude for domain shift'
    ENDIF
!---
    lon_temp(1:iim_file-i_cut-1) = lon_file(i_cut:iim_file)
    lon_temp(iim_file-i_cut:iim_file) = lon_file(1:i_cut+1)+360.
    lon_file(:) = lon_temp(:)
  ELSEIF (lon(1,1) < lon_file(1)) THEN
!-- A l'Ouest de la region d'interet, il n'y a pas de donnees
!-- le bout a l'Est de la region d'interet est deplace de 360 deg
!-- vers l'Ouest
    x_cut = lon(iim,1)
    DO i=1,iim_file
      IF (lon_file(i) < x_cut) i_cut = i
    ENDDO
    IF ( ( i_cut < 1 ) .OR. ( i_cut >= iim_file ) ) THEN
        STOP 'Cannot find longitude for domain shift'
    ENDIF
!---
    lon_temp(1:iim_file-i_cut-1) = lon_file(i_cut:iim_file) -360.
    lon_temp(iim_file-i_cut:iim_file) = lon_file(1:i_cut+1)
    lon_file(:) = lon_temp(:)
  ENDIF
!-
  DEALLOCATE (lon_temp)
  DEALLOCATE (lat_temp)
!-
  IF (    (lon(1,1) < lon_file(1)) &
 &    .OR.(     (lon(iim,1) > lon_file(iim_file)) &
 &         .AND.(lon(iim,1) > lon_file(1)+360.) ) ) THEN
    WRITE(numout,*) lon(:,1)
    WRITE(numout,*)
    WRITE(numout,*) lon_file(:)
    STOP 'weathgen_init: cannot find coherent x-coordinates'
  ENDIF
!-
  IF (i_cut /= 0) THEN
    CALL shift_field (iim_file,jjm_file,i_cut,lsm_file)
  ENDIF
!-
! 4. Verification
!-
  WRITE(numout,*)
  WRITE(numout,*) 'Input File: (Shifted) Global Land-Sea Mask'
  WRITE(numout,*)
  maskchar(0) = '-'
  maskchar(1) = 'X'
  plotstep = INT(REAL(iim_file-1)/termcol)+1
  DO j=1,jjm_file,plotstep
    DO i=1,iim_file,plotstep
      WRITE(numout,'(a1,$)') maskchar(NINT(lsm_file(i,j)))
    ENDDO
    WRITE(numout,*)
  ENDDO
  WRITE(numout,*)
!-
! 5. Prepare interpolation from fine grid land points to model grid
!-
! 5.1 how many grid points of the original grid fall into one grid
!     box of the model grid?
!-
  n_agg = NINT((zonal_res/xres_file*merid_res/yres_file )+1.)
!-
! au diable l'avarice !
!-
  n_agg = n_agg*2
!-
! 5.2 Allocate arrays where we store information about which
!     (and how many) points on the original grid fall
!     into which box on the model grid
!-

  ALLOC_ERR=-1
  ALLOCATE(ncorr(iim,jjm), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of ncorr : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(icorr(iim,jjm,n_agg), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of icorr : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(jcorr(iim,jjm,n_agg), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of jcorr : ",ALLOC_ERR
    STOP 
  ENDIF
!-
! 6. Land-Ocean Mask on model grid
!-
  tmp_lon = lon(:,1)
  tmp_lat = lat(1,:)

  CALL mask_c_o &
 &  (iim_file, jjm_file, lon_file, lat_file, lsm_file, fcrit, &
 &   iim, jjm, zonal_res, merid_res, n_agg, tmp_lon, tmp_lat, &
! &   iim, jjm, zonal_res, merid_res, n_agg, lon(:,1), lat(1,:), &
 &   mask, ncorr, icorr, jcorr)
!-
  WRITE(numout,*) 'Land-Sea Mask on Model Grid'
  WRITE(numout,*)
  plotstep = INT(REAL(iim-1)/termcol)+1
  DO j=1,jjm,plotstep
    DO i=1,iim,plotstep
      WRITE(numout,'(a1,$)') maskchar(NINT(mask(i,j)))
    ENDDO
    WRITE(numout,*)
  ENDDO
  WRITE(numout,*)
!-
! 7. kindex table.
!-
  nbindex = 0
  DO j=1,jjm
    DO i=1,iim
      IF (NINT(mask(i,j)) == 1) THEN
        nbindex = nbindex+1
        kindex(nbindex) = (j-1)*iim+i
      ENDIF
    ENDDO
  ENDDO
  nbindex_w = nbindex
  ALLOC_ERR=-1
  ALLOCATE(kindex_w(nbindex_w), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of kindex_w : ",ALLOC_ERR
    STOP 
  ENDIF
  kindex_w(:)=kindex(1:nbindex)
!-
  IF ( nbindex == 0 ) THEN
    WRITE(numout,*) 'Couillon! On est au plein milieu de l''ocean.'
    STOP 'Ou est-ce un bug?'
  ELSE
    WRITE(numout,*) 'Number of continental points: ',nbindex
  ENDIF

!-
! 10. clean up
!-
  DEALLOCATE (lon_file)
  DEALLOCATE (lat_file)
  DEALLOCATE (lsm_file)

END SUBROUTINE weathgen_init

!! ================================================================================================================================
!! SUBROUTINE   : weathgen_read_file
!!
!>\BRIEF         This subroutine allocates variables and reads files.  
!!
!! DESCRIPTION  :  
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): 
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE weathgen_read_file &
     & (force_id,iim,jjm)
  !---------------------------------------------------------------------
  IMPLICIT NONE
  !-

  !! 0. Variables and parameters declaration

  REAL,PARAMETER  :: fcrit = .5

  !! 0.1 Input variables 

  INTEGER,INTENT(IN)                  :: force_id
  INTEGER,INTENT(IN)                  :: iim, jjm

  !! 0.4 Local variables

  CHARACTER(LEN=30)               :: var_name

  INTEGER,DIMENSION(iim*jjm)  :: kindex
  INTEGER                     :: nbindex

  REAL,DIMENSION(iim*jjm)     :: xchamp
  REAL,DIMENSION(iim*jjm,nmon)  :: xchampm

!_ ================================================================================================================================

  kindex(:)=0

#ifdef CPP_PARA
  nbindex=nbp_loc
  CALL scatter(index_g,kindex)
  kindex(1:nbindex)=kindex(1:nbindex)-(jj_begin-1)*iim_g
#else
  ! Copy saved land points index
  nbindex = nbindex_w
  kindex(1:nbindex_w) = kindex_w(:)
#endif

!-

  ALLOC_ERR=-1
  ALLOCATE(lat_land(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of lat_land : ",ALLOC_ERR
    STOP 
  ENDIF

!-
! 8 topography
!-

  ALLOC_ERR=-1
  ALLOCATE(xintopo(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of xintopo : ",ALLOC_ERR
    STOP 
  ENDIF
!-
  var_name='altitude'
  CALL weather_read (force_id,var_name,iim_file,jjm_file,1,i_cut, &
                     iim,jjm,n_agg,ncorr,icorr,jcorr,xchamp)
  xintopo(:)=xchamp(kindex(1:nbindex))
!-
! 9. Allocate and read the monthly fields
!-

  ALLOC_ERR=-1
  ALLOCATE(xinwet(nbindex,nmon), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of xinwet : ",ALLOC_ERR
    STOP 
  ENDIF
  var_name='prs'
  CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
                     iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
  xinwet(:,:)=xchampm(kindex(1:nbindex),:)
  !-

  ALLOC_ERR=-1
  ALLOCATE(xinprec(nbindex,nmon), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of xinprec : ",ALLOC_ERR
    STOP 
  ENDIF
  var_name='prm'
  CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
                     iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
  xinprec(:,:)=xchampm(kindex(1:nbindex),:)
!-
  
  ALLOC_ERR=-1
  ALLOCATE(xint(nbindex,nmon), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of xint : ",ALLOC_ERR
    STOP 
  ENDIF
  var_name='t2m'
  CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
                     iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
  xint(:,:)=xchampm(kindex(1:nbindex),:)
!-

  ALLOC_ERR=-1
  ALLOCATE(xinq(nbindex,nmon), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of xinq : ",ALLOC_ERR
    STOP 
  ENDIF
  var_name='r2m'
  CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
                     iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
  xinq(:,:)=xchampm(kindex(1:nbindex),:)
!-

  ALLOC_ERR=-1
  ALLOCATE(xinwind(nbindex,nmon), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of xinwind : ",ALLOC_ERR
    STOP 
  ENDIF
  var_name='uv10m'
  CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
                     iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
  xinwind(:,:)=xchampm(kindex(1:nbindex),:)
!-

  ALLOC_ERR=-1
  ALLOCATE(xincld(nbindex,nmon), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of xincld : ",ALLOC_ERR
    STOP 
  ENDIF
  var_name='tc'
  CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
                       iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
  xincld(:,:)=xchampm(kindex(1:nbindex),:)
!-

  ALLOC_ERR=-1
  ALLOCATE(xintrng(nbindex,nmon), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of xintrng : ",ALLOC_ERR
    STOP 
  ENDIF
  var_name='t2a'
  CALL weather_read (force_id,var_name,iim_file,jjm_file,nmon,i_cut, &
                       iim,jjm,n_agg,ncorr,icorr,jcorr,xchampm)
  xintrng(:,:)=xchampm(kindex(1:nbindex),:)
!-
! 10. clean up
!-
  IF (is_root_prc) THEN
     DEALLOCATE (ncorr)
     DEALLOCATE (icorr)
     DEALLOCATE (jcorr)
  ENDIF

!-
! 12. Allocate space for daily mean values
!-

  ALLOC_ERR=-1
  ALLOCATE(iwet(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of iwet : ",ALLOC_ERR
    STOP 
  ENDIF
!-

  ALLOC_ERR=-1
  ALLOCATE(psurfm0(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of psurfm0 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(cloudm0(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of cloudm0 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(tmaxm0(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of tmaxm0 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(tminm0(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of tminm0 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(qdm0(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of qdm0 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(udm0(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of udm0 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(precipm0(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of precipm0 : ",ALLOC_ERR
    STOP 
  ENDIF
!-

  ALLOC_ERR=-1
  ALLOCATE(psurfm1(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of psurfm1 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(cloudm1(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of cloudm1 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(tmaxm1(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of tmaxm1 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(tminm1(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of tminm1 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(qdm1(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of qdm1 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(udm1(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of udm1 : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(precipm1(nbindex), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of precipm1 : ",ALLOC_ERR
    STOP 
  ENDIF
END SUBROUTINE weathgen_read_file
!
!=
!

!! ================================================================================================================================
!! SUBROUTINE   : weathgen_begin
!!
!>\BRIEF         This subroutines initializes variables or reads restart values. 
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::nbindex, ::tair, ::pb, ::qair, ::kindex
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE weathgen_begin ( &
     & dt_force,itau, date0, &
     & rest_id,iim,jjm, &
     & lon,lat,tair,pb,qair,kindex,nbindex)
  !---------------------------------------------------------------------
  IMPLICIT NONE

  !-

  !! 0. Variables and parameters declaration

  !! 0.1 Input variables

  REAL,INTENT(IN)                     :: dt_force, date0
  INTEGER,INTENT(IN)                  :: itau
  INTEGER,INTENT(IN)                  :: rest_id
  INTEGER,INTENT(IN)                  :: iim, jjm
  REAL,DIMENSION(iim,jjm),INTENT(IN)  :: lon,lat

  !! 0.2 Output variables 

  INTEGER,INTENT(OUT)                     :: nbindex
  INTEGER,DIMENSION(iim*jjm),INTENT(OUT)  :: kindex
  REAL,DIMENSION(iim,jjm),INTENT(OUT)     :: tair,pb,qair

  !! 0.4 Local variables

  INTEGER  :: i, j, ij
  REAL     :: val_exp
  REAL,DIMENSION(iim*jjm)               :: xchamp
  REAL,DIMENSION(iim_g*jjm_g)           :: xchamp_g
  CHARACTER(LEN=30)                     :: var_name
  REAL,DIMENSION(1)                     :: jullasttab
  REAL,DIMENSION(seedsize_max)          :: seed_in_file
  INTEGER,DIMENSION(:), ALLOCATABLE     :: seed_in_proc
  INTEGER  :: seedsize, iret
  INTEGER  :: nlonid, nlatid, nlonid1, nlatid1, tdimid1, varid
  INTEGER  :: ndim, dims(4)
  CHARACTER(LEN=30) :: assoc
  CHARACTER(LEN=70) :: str70
  CHARACTER(LEN=80) :: stamp
  INTEGER  :: yy_b, mm_b, dd_b, hh, mm
  REAL     :: ss
  CHARACTER(LEN=10) :: today, att
  INTEGER  :: nlandid1, nlandid, nlevid, nlevid1
  REAL     :: lev_max, lev_min
  REAL     :: height_lev1
  INTEGER  :: imois
  REAL     :: xx, td

!_ ================================================================================================================================

  kindex(:)=0

#ifdef CPP_PARA
  nbindex=nbp_loc
  CALL scatter(index_g,kindex)
  kindex(1:nbindex)=kindex(1:nbindex)-(jj_begin-1)*iim_g
#else
  ! Copy saved land points index
  nbindex = nbindex_w
  kindex(1:nbindex_w) = kindex_w(:)
#endif
!-
! 13. Prescribed or statistical values?
!-
  !Config  Key  = IPPREC
  !Config  Desc = Use prescribed values
  !Config  If = ALLOW_WEATHERGEN
  !Config  Def  = 0
  !Config  Help = If this is set to 1, the weather generator
  !Config         uses the monthly mean values for daily means.
  !Config         If it is set to 0, the weather generator
  !Config         uses statistical relationships to derive daily
  !Config         values from monthly means.
  ipprec = 0
  CALL getin_p ('IPPREC',ipprec)
  WRITE(numout,*) 'IPPREC: ',ipprec
!-
! 14. Do we want exact monthly precipitations even with ipprec=0 ?
!-
  !Config  Key  = WEATHGEN_PRECIP_EXACT
  !Config  Desc = Exact monthly precipitation
  !Config  If = ALLOW_WEATHERGEN
  !Config  Def  = n
  !Config  Help = If this is set to y, the weather generator
  !Config         will generate pseudo-random precipitations
  !Config         whose monthly mean is exactly the prescribed one.
  !Config         In this case, the daily precipitation (for rainy
  !Config         days) is constant (that is, some days have 0 precip,
  !Config         the other days have precip=Precip_month/n_precip,
  !Config         where n_precip is the prescribed number of rainy days
  !Config         per month).
  precip_exact = .FALSE.
  CALL getin_p ('WEATHGEN_PRECIP_EXACT',precip_exact)
  WRITE(numout,*) 'PRECIP_EXACT: ',precip_exact
!-
  IF (precip_exact) THEN
!---
!-- preparer un tableau utilise pour determiner s'il pleut ou pas
!-- (en fct. du nombre de jours de pluie par mois).
!---
     IF (is_root_prc) THEN
        DO imois=1,12
           CALL permute (ndaypm(imois),jour_precip(:,imois))
        ENDDO
     ENDIF
     CALL bcast(jour_precip)
  ENDIF
!-
! Read Orbital Parameters
!-
  !Config  Key  = ECCENTRICITY
  !Config  Desc = Use prescribed values
  !Config  If = ALLOW_WEATHERGEN
  !Config  Def  = 0.016724
  ecc = 0.016724
  CALL getin_p ('ECCENTRICITY',ecc)
  WRITE(numout,*) 'ECCENTRICITY: ',ecc
!
  !Config  Key  = PERIHELIE
  !Config  Desc = Use prescribed values
  !Config  If = ALLOW_WEATHERGEN
  !Config  Def  = 102.04
  perh = 102.04
  CALL getin_p ('PERIHELIE',perh)
  WRITE(numout,*) 'PERIHELIE: ',perh
!
  !Config  Key  = OBLIQUITY
  !Config  Desc = Use prescribed values
  !Config  If = ALLOW_WEATHERGEN
  !Config  Def  = 23.446
  xob = 23.446
  CALL getin_p ('OBLIQUITY',xob)
  WRITE(numout,*) 'OBLIQUITY: ',xob
!-
! 15. Read restart file
!-
  CALL ioget_expval (val_exp)
!-
  var_name= 'julian'
  IF (is_root_prc) THEN
     CALL restget (rest_id,var_name,1,1,1,itau,.TRUE.,jullasttab)
     IF (jullasttab(1) == val_exp) THEN
        jullasttab(1) = itau2date(itau-1, date0, dt_force)
     ENDIF
  ENDIF
  CALL bcast(jullasttab)
  julian_last = jullasttab(1)
!-
  var_name= 'seed'
  IF (is_root_prc) &
       CALL restget (rest_id,var_name,seedsize_max, &
 &              1,1,itau,.TRUE.,seed_in_file)
  CALL bcast(seed_in_file)
  IF (ALL(seed_in_file(:) == val_exp)) THEN
!---
!-- there is no need to reinitialize the random number generator as
!-- this does not seem to be a restart
!---
    CONTINUE
  ELSE
!---
!-- reinitialize the random number generator
!---
     IF (is_root_prc) &
          CALL RANDOM_SEED( SIZE = seedsize )
     CALL bcast(seedsize)
    IF (seedsize > seedsize_max) THEN
        STOP 'weathgen_begin: increase seedsize_max'
    ENDIF
  
    ALLOC_ERR=-1
    ALLOCATE(seed_in_proc(seedsize), STAT=ALLOC_ERR)
    IF (ALLOC_ERR/=0) THEN
      WRITE(numout,*) "ERROR IN ALLOCATION of seed_in_proc : ",ALLOC_ERR
      STOP 
    ENDIF
    seed_in_proc(1:seedsize) = NINT( seed_in_file(1:seedsize) )
    CALL RANDOM_SEED (PUT = seed_in_proc)
    DEALLOCATE( seed_in_proc )
  ENDIF
!-
  var_name= 'iwet'
  IF (is_root_prc) THEN
    CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
    IF (ALL(xchamp_g(:) == val_exp)) THEN
      xchamp_g(:) = 0.
    ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  iwet(:) = NINT(xchamp(kindex(1:nbindex)))
!-
  var_name= 'psurfm0'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 100000.
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  psurfm0(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'cloudm0'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = .5
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  cloudm0(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'tmaxm0'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 285.
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  tmaxm0(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'tminm0'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 275.
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  tminm0(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'qdm0'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 1.E-03
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  qdm0(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'udm0'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 2.
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  udm0(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'precipm0'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 1.
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  precipm0(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'psurfm1'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 100000.
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  psurfm1(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'cloudm1'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = .5
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  cloudm1(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'tmaxm1'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 285.
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  tmaxm1(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'tminm1'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 275.
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  tminm1(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'qdm1'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 1.E-03
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  qdm1(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'udm1'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 2.
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  udm1(:) = xchamp(kindex(1:nbindex))
!-
  var_name= 'precipm1'
  IF (is_root_prc) THEN
     CALL restget (rest_id, var_name, iim_g, jjm_g, 1, itau, .TRUE., xchamp_g)
     IF (ALL(xchamp_g(:) == val_exp)) THEN
        xchamp_g(:) = 1.
     ENDIF
  ENDIF
  CALL scatter2D(xchamp_g,xchamp)
  precipm1(:) = xchamp(kindex(1:nbindex))
!-
! 16. We still need instantaneous tair, qair, and the surface pressure
!     We take daily mean values read from the restart file
!-
!!$  tair(:,:)=280.
!!$  qair(:,:)=1.E-03
!!$  pb(:,:)=101325
  tair(:,:)=val_exp
  qair(:,:)=val_exp
  pb(:,:)=val_exp
  xx = 9.81/287./0.0065
  DO ij=1,nbindex
     j = ((kindex(ij)-1)/iim) + 1
     i = kindex(ij) - (j-1)*iim
     
     lat_land(ij) = lat(i,j)
     
     td =  (tmaxm0(ij)+tminm0(ij))/2.
     tair(i,j) = td
     qair(i,j) = qdm1(ij)
     pb(i,j) = 101325.*(td/(td+0.0065*xintopo(ij)))**xx
  ENDDO
!-
! 17. We can write a forcing file for Orchidee
!     from this weather Generator run.
!-
  !Config  Key  = DUMP_WEATHER
  !Config  Desc = Write weather from generator into a forcing file
  !Config  Def  = n
  !Config  Help = This flag makes the weather generator dump its
  !               generated weather into a forcing file which can
  !               then be used to get the same forcing on different
  !               machines. This only works correctly if there is
  !               a restart file (otherwise the forcing at the first
  !               time step is slightly wrong).
  dump_weather = .FALSE.
  CALL getin_p ('DUMP_WEATHER',dump_weather)
!-
  IF (dump_weather .AND. is_root_prc) THEN
!---
!-- Initialize the file
!---
    !Config  Key  = DUMP_WEATHER_FILE
    !Config  Desc = Name of the file that contains
    !               the weather from generator
    !Config  Def  = 'weather_dump.nc'
    !Config  If = DUMP_WEATHER
    !Config  Help =
    dump_weather_file = 'weather_dump.nc'
    CALL getin ('DUMP_WEATHER_FILE',dump_weather_file)
!---
    !Config  Key  = DUMP_WEATHER_GATHERED
    !Config  Desc = Dump weather data on gathered grid
    !Config  Def  = y
    !Config  If = DUMP_WEATHER
    !Config  Help = If 'y', the weather data are gathered
    !               for all land points.
    gathered = .TRUE.
    CALL getin ('DUMP_WEATHER_GATHERED',gathered)
!---
    iret = NF90_CREATE (TRIM(dump_weather_file),NF90_CLOBBER,dump_id)
!---
!-- Dimensions
!---
    iret = NF90_DEF_DIM (dump_id,'x',iim_g,nlonid1)
    iret = NF90_DEF_DIM (dump_id,'y',jjm_g,nlatid1)
    iret = NF90_DEF_DIM (dump_id,'z',  1,nlevid1)
!---
    IF (gathered) THEN
      iret = NF90_DEF_DIM (dump_id,'land',nbp_glo,nlandid1)
    ENDIF
    iret = NF90_DEF_DIM (dump_id,'tstep',NF90_UNLIMITED,tdimid1)
!---
!-- Coordinate variables
!---
    dims(1:2) = (/ nlonid1, nlatid1 /)
!---
    iret = NF90_DEF_VAR (dump_id,'nav_lon',n_rtp,dims(1:2),nlonid)
    iret = NF90_PUT_ATT (dump_id,nlonid,'units',"degrees_east")
    iret = NF90_PUT_ATT (dump_id,nlonid,'valid_min',MINVAL(lon_g))
    iret = NF90_PUT_ATT (dump_id,nlonid,'valid_max',MAXVAL(lon_g))
    iret = NF90_PUT_ATT (dump_id,nlonid,'long_name',"Longitude")
!---
    iret = NF90_DEF_VAR (dump_id,'nav_lat',n_rtp,dims(1:2),nlatid)
    iret = NF90_PUT_ATT (dump_id,nlatid,'units',"degrees_north")
    iret = NF90_PUT_ATT (dump_id,nlatid,'valid_min',MINVAL(lat_g))
    iret = NF90_PUT_ATT (dump_id,nlatid,'valid_max',MAXVAL(lat_g))
    iret = NF90_PUT_ATT (dump_id,nlatid,'long_name',"Latitude")
!---
    height_lev1 = 10.
    CALL getin ('HEIGHT_LEV1',height_lev1)
    lev_min = height_lev1
    lev_max = height_lev1
!---
    iret = NF90_DEF_VAR (dump_id,'level',n_rtp,(/ nlevid1 /),nlevid)
    iret = NF90_PUT_ATT (dump_id,nlevid,'units',"m")
    iret = NF90_PUT_ATT (dump_id,nlevid,'valid_min',lev_min)
    iret = NF90_PUT_ATT (dump_id,nlevid,'valid_max',lev_max)
    iret = NF90_PUT_ATT (dump_id,nlevid,'long_name',"Vertical levels")
!---
    IF (gathered) THEN
      iret = NF90_DEF_VAR (dump_id,'land',NF90_INT,(/ nlandid1 /),nlandid)
      iret = NF90_PUT_ATT (dump_id,nlandid,'compress',"y x")
    ENDIF
!---
!-- Store the time axes.
!---
    iret = NF90_DEF_VAR (dump_id,'time',n_rtp,tdimid1,time_id)

    yy_b=0
    mm_b=1
    dd_b=1
    hh=00
    mm=00
    ss=0.

    WRITE (str70,7000) yy_b, mm_b, dd_b, hh, mm, INT(ss)
    iret = NF90_PUT_ATT (dump_id,time_id,'units',TRIM(str70))
    iret = NF90_PUT_ATT (dump_id,time_id,'calendar',TRIM(calendar_str))
    iret = NF90_PUT_ATT (dump_id,time_id,'title','Time')
    iret = NF90_PUT_ATT (dump_id,time_id,'long_name','Time axis')
    WRITE(str70,7001) yy_b, cal(mm_b), dd_b, hh, mm, INT(ss)
    iret = NF90_PUT_ATT (dump_id,time_id,'time_origin',TRIM(str70))
!---
!-- Time steps
!---
    iret = NF90_DEF_VAR (dump_id,'timestp',NF90_INT,tdimid1,timestp_id)
    WRITE(str70,7002) yy_b, mm_b, dd_b, hh, mm, INT(ss)
    iret = NF90_PUT_ATT (dump_id,timestp_id,'units',TRIM(str70))
    iret = NF90_PUT_ATT (dump_id,timestp_id,'title','Time steps')
    iret = NF90_PUT_ATT (dump_id,timestp_id,'tstep_sec',dt_force)
    iret = NF90_PUT_ATT &
 &           (dump_id,timestp_id,'long_name','Time step axis')
    WRITE(str70,7001) yy_b, cal(mm_b), dd_b, hh, mm, INT(ss)
    iret = NF90_PUT_ATT (dump_id,timestp_id,'time_origin',TRIM(str70))
!---
7000 FORMAT('seconds since ',I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)
7001 FORMAT(' ',I4.4,'-',A3,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)
7002 FORMAT('timesteps since ', I4.4,'-',I2.2,'-',I2.2,' ',I2.2,':',I2.2,':',I2.2)
!---
!-- The variables in the file are declared and their attributes
!-- written.
!---
    IF (gathered) THEN
      ndim = 2
      dims(1:2) = (/ nlandid1, tdimid1 /)
      assoc = 'time (nav_lat nav_lon)'
    ELSE
      ndim = 3
      dims(1:3) = (/ nlonid1, nlatid1, tdimid1 /)
      assoc = 'time nav_lat nav_lon'
    ENDIF
!---
    iret = NF90_DEF_VAR (dump_id,'SWdown',n_rtp,dims(1:ndim),varid)
    iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX')
    iret = NF90_PUT_ATT (dump_id,varid,'units','W/m^2')
    iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
 &                       'Surface incident shortwave radiation')
    iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
    iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef)
    soldownid = varid
!---
    iret = NF90_DEF_VAR (dump_id,'Rainf',n_rtp,dims(1:ndim),varid)
    iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX')
    iret = NF90_PUT_ATT (dump_id,varid,'units','kg/m^2s')
    iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
 &                       'Rainfall rate')
    iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
    iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef)
    rainfid = varid
!---
    iret = NF90_DEF_VAR (dump_id,'Snowf',n_rtp,dims(1:ndim),varid)
    iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX')
    iret = NF90_PUT_ATT (dump_id,varid,'units','kg/m^2s')
    iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
 &                       'Snowfall rate')
    iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
    iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef)
    snowfid = varid
!---
    iret = NF90_DEF_VAR (dump_id,'LWdown',n_rtp,dims(1:ndim),varid)
    iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX')
    iret = NF90_PUT_ATT (dump_id,varid,'units','W/m^2')
    iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
 &                       'Surface incident longwave radiation')
    iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
    iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef)
    lwradid = varid
!---
    iret = NF90_DEF_VAR (dump_id,'PSurf',n_rtp,dims(1:ndim),varid)
    iret = NF90_PUT_ATT (dump_id,varid,'axis','TYX')
    iret = NF90_PUT_ATT (dump_id,varid,'units','Pa')
    iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
 &                       'Surface pressure')
    iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
    iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef)
    psolid = varid
!---
!-- 3D Variables to be written
!---
    IF (gathered) THEN
      ndim = 3
      dims(1:3) = (/ nlandid1, nlevid1, tdimid1 /)
      assoc = 'time level (nav_lat nav_lon)'
    ELSE
      ndim = 4
      dims(1:4) = (/ nlonid1, nlatid1, nlevid1, tdimid1 /)
      assoc = 'time level nav_lat nav_lon'
    ENDIF
!---
    iret = NF90_DEF_VAR (dump_id,'Tair',n_rtp,dims(1:ndim),varid)
    iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX')
    iret = NF90_PUT_ATT (dump_id,varid,'units','K')
    iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
 &                       'Near surface air temperature')
    iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
    iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef)
    tairid = varid
!---
    iret = NF90_DEF_VAR (dump_id,'Qair',n_rtp,dims(1:ndim),varid)
    iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX')
    iret = NF90_PUT_ATT (dump_id,varid,'units','kg/kg')
    iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
 &                       'Near surface specific humidity')
    iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
    iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef)
    qairid = varid
!---
    iret = NF90_DEF_VAR (dump_id,'Wind_N',n_rtp,dims(1:ndim),varid)
    iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX')
    iret = NF90_PUT_ATT (dump_id,varid,'units','m/s')
    iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
 &                       'Near surface northward wind component')
    iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
    iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef)
    uid = varid
!---
    iret = NF90_DEF_VAR (dump_id,'Wind_E',n_rtp,dims(1:ndim),varid)
    iret = NF90_PUT_ATT (dump_id,varid,'axis','TZYX')
    iret = NF90_PUT_ATT (dump_id,varid,'units','m/s')
    iret = NF90_PUT_ATT (dump_id,varid,'long_name', &
 &                       'Near surface eastward wind component')
    iret = NF90_PUT_ATT (dump_id,varid,'associate',TRIM(assoc))
    iret = NF90_PUT_ATT (dump_id,varid,'missing_value',undef)
    vid = varid
!---
!-- Global attributes
!---
    CALL DATE_AND_TIME (today, att)
    stamp = "WG, date: "//TRIM(today)//" at "//TRIM(att)
    iret = NF90_PUT_ATT (dump_id,NF90_GLOBAL,'Conventions',"GDT 1.2")
    iret = NF90_PUT_ATT (dump_id,NF90_GLOBAL,'file_name', &
 &                       TRIM(dump_weather_file))
    iret = NF90_PUT_ATT (dump_id,NF90_GLOBAL,'production',TRIM(stamp))
!---
!-- Finish the definition phase
!---
    iret = NF90_ENDDEF (dump_id)
!---
!--    Write coordinates
!---
    iret = NF90_PUT_VAR (dump_id,nlonid,lon_g)
    IF (iret /= NF90_NOERR) THEN
       WRITE(numout,*) iret
       CALL ipslerr (3,'weathgen_begin', &
 &          'Could not put variable nav_lon in the file : ', &
 &          TRIM(dump_weather_file),'(Solution ?)')
    ENDIF
    iret = NF90_PUT_VAR (dump_id,nlatid,lat_g)
    IF (iret /= NF90_NOERR) THEN
       WRITE(numout,*) iret
       CALL ipslerr (3,'weathgen_begin', &
 &          'Could not put variable nav_lat in the file : ', &
 &          TRIM(dump_weather_file),'(Solution ?)')
    ENDIF
    iret = NF90_PUT_VAR (dump_id,nlevid,height_lev1)
    IF (iret /= NF90_NOERR) THEN
       WRITE(numout,*) iret
       CALL ipslerr (3,'weathgen_begin', &
 &          'Could not put variable level in the file : ', &
 &          TRIM(dump_weather_file),'(Solution ?)')
    ENDIF
!---
    IF (gathered) THEN
       iret = NF90_PUT_VAR (dump_id,nlandid,index_g(1:nbp_glo))
       IF (iret /= NF90_NOERR) THEN
          WRITE(numout,*) iret
          CALL ipslerr (3,'weathgen_begin', &
 &          'Could not put variable land in the file : ', &
 &          TRIM(dump_weather_file),'(Solution ?)')
       ENDIF
    ENDIF
!---
 ENDIF ! dump_weather
!-----------------------------
END SUBROUTINE weathgen_begin
!-
!===
!-
!! ================================================================================================================================
!! SUBROUTINE   : weathgen_get
!!
!>\BRIEF        This subroutine initializes forcing values. 
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::swdown, ::raina, ::snowa, ::tair, ::u, ::v, ::qair, ::psurf, ::lwdown
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE weathgen_get &
& (itau, date0, dt_force, nbindex, nband, lat, &
&  swdown, raina, snowa, tair, u, v, qair, psurf, lwdown)
!---------------------------------------------------------------------
  IMPLICIT NONE

  !! 0. Variables and parameters declaration

  !! 0.1 Input variables

  INTEGER,INTENT(IN)               :: itau      !! number of time step
  REAL,INTENT(IN)                  :: date0     !! date when itau was 0
  REAL,INTENT(IN)                  :: dt_force  !! time step (s)
  INTEGER,INTENT(IN)               :: nbindex   !! number of land points
  INTEGER,INTENT(IN)               :: nband     !! number of visible bands
  REAL,DIMENSION(nbindex),INTENT(IN)  :: lat    !! latitude (deg)

  !! 0.2 Output variables  

  REAL,DIMENSION(nbindex,nband),INTENT(OUT) :: swdown        !! Downward solar radiation @tex $(W.m^{-2})$ @endtex
  REAL,DIMENSION(nbindex),INTENT(OUT)       :: raina, snowa  !! precipitations @tex $(mm.day^{-1})$ @endtex 
  REAL,DIMENSION(nbindex),INTENT(OUT)       :: tair          !! air temperature (K)
  REAL,DIMENSION(nbindex),INTENT(OUT)       :: u,v           !! wind speed @tex $(m.s^{-1})$ @endtex
  REAL,DIMENSION(nbindex),INTENT(OUT)       :: qair          !! air humidity @tex $(kg_h2o.kg_air^{-1})$ @endtex
  REAL,DIMENSION(nbindex),INTENT(OUT)       :: psurf         !! Surface pressure (Pa)
  REAL,DIMENSION(nbindex),INTENT(OUT)       :: lwdown        !! Downward long wave radiation @tex $(W.m^{-2})$ @endtex

  !! 0.4 Local variables

  REAL,DIMENSION(nbindex)        :: cloud, tmax, tmin, precipd, qd, ud
  REAL,DIMENSION(nbindex)        :: rh
  REAL,DIMENSION(nbindex,nband)  :: solai, solad
  REAL                        :: julian, jur
  REAL                        :: x
  INTEGER                     :: yy, mm, dd
  REAL                        :: ss, plens, time

!_ ================================================================================================================================

!-
! 1. get a reduced julian day
!-
  julian = itau2date(itau-1, date0, dt_force)
!SZ, test: solar noon at 12 o'clock!
!  julian = itau2date(itau, date0, dt_force)
  CALL ju2ymds (julian, yy, mm, dd, ss)
  CALL ymds2ju (yy,1,1,0.0, jur)
  julian = julian-jur
  CALL ju2ymds (julian, yy, mm, dd, ss)
!-
! 2. daily values
!-
  IF (INT(julian_last) /= INT(julian)) THEN
!-- today's values become yesterday's values
    cloudm1(:) = cloudm0(:)
    tmaxm1(:) = tmaxm0(:)
    tminm1(:) = tminm0(:)
    precipm1(:) = precipm0(:)
    qdm1(:) = qdm0(:)
    udm1(:) = udm0(:)
    psurfm1(:) = psurfm0(:)
!-- we have to get new daily values
!!$    WRITE(*,*) mpi_rank, "weathgen_get : date ",yy, mm, dd, ss
!!$    WRITE(*,*) mpi_rank, "weathgen_get : grid date ",year, month, day, sec
    CALL daily (nbindex, mm, dd, cloudm0, tmaxm0, tminm0, &
 &              precipm0, qdm0, udm0, psurfm0)
  ENDIF
!-
! 3. interpolate daily values
!    (otherwise we get ugly temperature jumps at midnight)
!-
  x = (julian-INT(julian))
!-
  cloud(:) = (1.-x)*cloudm1(:)+x*cloudm0(:)
  tmax(:) = (1.-x)*tmaxm1(:)+x*tmaxm0(:)
  tmin(:) = (1.-x)*tminm1(:)+x*tminm0(:)
  precipd(:) = (1.-x)*precipm1(:)+x*precipm0(:)
  qd(:) = (1.-x)*qdm1(:)+x*qdm0(:)
  ud(:) = (1.-x)*udm1(:)+x*udm0(:)
  psurf(:) = (1.-x)*psurfm1(:)+x*psurfm0(:)
!-
! 4. read instantaneous values
!-
  plens = one_day/dt_force
  time = (julian-REAL(INT(julian)))*one_day
!-
  CALL diurnal (nbindex, nband, time, NINT(julian), plens, 0., one_day, &
 &              lat, cloud, tmax, tmin, precipd, qd, ud, psurf, &
 &              lwdown, solad, solai, u, tair, qair, raina, snowa, rh)
!-
  raina(:) = raina(:)/dt_force
  snowa(:) = snowa(:)/dt_force
!-
  swdown(:,:) = solad(:,:)+solai(:,:)
!-
  v(:) = 0.
!-
! 5. Store date
!-
  julian_last = julian
!--------------------------
END SUBROUTINE weathgen_get
!-
!===
!-
!! ================================================================================================================================
!! SUBROUTINE   : weathgen_restwrite
!!
!>\BRIEF         This subroutine writes data in the restart file.
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): None
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE weathgen_restwrite (rest_id,itau,iim,jjm,nbindex,kindex)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-

  !! 0. Variables and parameters declaration

  !! 0.1 Input variables

  INTEGER,INTENT(IN)  :: rest_id,itau,iim,jjm,nbindex
  INTEGER,DIMENSION(iim*jjm),INTENT(IN) :: kindex

  !! 0.4 Local variables

  CHARACTER(LEN=30)                 :: var_name
  INTEGER                           :: i,j,ij
  REAL,DIMENSION(1)                 :: jullasttab
  REAL,DIMENSION(seedsize_max)      :: seed_in_file
  INTEGER,DIMENSION(:),ALLOCATABLE  :: seed_in_proc
  INTEGER                           :: seedsize
  REAL                              :: rndnum
  REAL,DIMENSION(iim*jjm)           :: xchamp
  REAL,DIMENSION(iim_g*jjm_g)       :: xchamp_g

!_ ================================================================================================================================

  var_name= 'julian'
  jullasttab(:) = julian_last
  IF (is_root_prc) CALL restput (rest_id, var_name, 1, 1, 1, itau, jullasttab)
!-
  IF (is_root_prc) THEN
     CALL RANDOM_SEED( SIZE = seedsize )
     IF (seedsize > seedsize_max) THEN
        STOP 'weathgen_restwrite: increase seedsize_max'
     ENDIF
  ENDIF
  CALL bcast(seedsize)

  IF (is_root_prc) THEN
     ALLOC_ERR=-1
     ALLOCATE(seed_in_proc(seedsize), STAT=ALLOC_ERR)
     IF (ALLOC_ERR/=0) THEN
        WRITE(numout,*) "ERROR IN ALLOCATION of seed_in_proc : ",ALLOC_ERR
        STOP 
     ENDIF
!-
     CALL RANDOM_SEED (GET = seed_in_proc)
!-
     seed_in_file(1:seedsize) = REAL(seed_in_proc(1:seedsize))
!-
! fill in the seed up to seedsize_max
! (useful in the case we restart on
!  a machine with a longer seed vector:
!  we do not want a degenerated seed)
!-
     DO i=seedsize+1,seedsize_max
        CALL RANDOM_NUMBER( rndnum )
        seed_in_file(i) = 100000.*rndnum
     ENDDO
  ENDIF
  CALL bcast (seed_in_file)
!-
  IF (is_root_prc) THEN
     DEALLOCATE( seed_in_proc )
!-
     var_name= 'seed'
     CALL restput (rest_id,var_name,seedsize_max,1,1,itau,seed_in_file)
  ENDIF
!-

  xchamp(:) = val_exp
  
!!$  DO j=1,jjm
!!$    DO i=1,iim
!!$      ij = (j-1)*iim+i
!!$      xchamp(i,j) = REAL(iwet(ij))
!!$    ENDDO
!!$  ENDDO
  DO ij=1,nbindex
     xchamp(kindex(ij)) = REAL(iwet(ij))
  ENDDO
  var_name= 'iwet'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = psurfm0(ij)
  ENDDO
  var_name= 'psurfm0'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = cloudm0(ij)
  ENDDO
  var_name= 'cloudm0'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = tmaxm0(ij)
  ENDDO
  var_name= 'tmaxm0'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = tminm0(ij)
  ENDDO
  var_name= 'tminm0'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = qdm0(ij)
  ENDDO
  var_name= 'qdm0'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = udm0(ij)
  ENDDO
  var_name= 'udm0'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = precipm0(ij)
  ENDDO
  var_name= 'precipm0'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = psurfm1(ij)
  ENDDO
  var_name= 'psurfm1'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = cloudm1(ij)
  ENDDO
  var_name= 'cloudm1'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = tmaxm1(ij)
  ENDDO
  var_name= 'tmaxm1'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = tminm1(ij)
  ENDDO
  var_name= 'tminm1'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = qdm1(ij)
  ENDDO
  var_name= 'qdm1'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = udm1(ij)
  ENDDO
  var_name= 'udm1'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!-
  DO ij=1,nbindex
    xchamp(kindex(ij)) = precipm1(ij)
  ENDDO
  var_name= 'precipm1'
  CALL gather2D(xchamp,xchamp_g)
  IF (is_root_prc) THEN
     CALL restput (rest_id, var_name, iim_g, jjm_g, 1, itau, xchamp_g)
  ENDIF
!--------------------------------
END SUBROUTINE weathgen_restwrite
!-
!===
!-
!! ================================================================================================================================
!! SUBROUTINE   : weathgen_data
!!
!>\BRIEF         This subroutine reads data from an netcdf file. 
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::champout
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE weather_read &
& (force_id,nomvar,iim_file,jjm_file,n3,i_cut, &
&  iim,jjm,n_agg,ncorr,icorr,jcorr,champout)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-

  !! 0. Variables and parameters declaration
  
  !! 0.1 Input variables

  INTEGER,INTENT(IN)                          :: force_id
  CHARACTER(LEN=*),INTENT(IN)                 :: nomvar
  INTEGER,INTENT(IN)                          :: iim_file,jjm_file
  INTEGER,INTENT(IN)                          :: n3
  INTEGER,INTENT(IN)                          :: i_cut
  INTEGER,INTENT(IN)                          :: iim,jjm
  INTEGER,INTENT(IN)                          :: n_agg
  INTEGER,DIMENSION(:,:),INTENT(IN)           :: ncorr
  INTEGER,DIMENSION(:,:,:),INTENT(IN)         :: icorr,jcorr

  !! 0.2 Output variables 

  REAL,DIMENSION(iim*jjm,n3),INTENT(OUT)      :: champout

  !! 0.4 Local variables

  REAL,DIMENSION(iim_file,jjm_file,n3)        :: champ_file
  REAL,ALLOCATABLE,DIMENSION(:,:)             :: champout_g
  INTEGER                                     :: i,j,ij,l,m

!_ ================================================================================================================================

  WRITE(numout,*) 'Lecture ',TRIM(nomvar)
!-
  IF (is_root_prc) THEN
     ALLOCATE(champout_g(iim_g*jjm_g,n3))
     IF ( n3 == 1 ) THEN
        CALL flinget (force_id,nomvar(1:LEN_TRIM(nomvar)), &
             &                iim_file, jjm_file, 0, 0,  1, 1, champ_file)
     ELSE
        DO l=1,n3
           CALL flinget &
                &      (force_id,nomvar(1:LEN_TRIM(nomvar)), &
                &       iim_file, jjm_file, 0, n3,  l, l, champ_file(:,:,l))
        ENDDO
     ENDIF
     ! shift if necessary
     IF (i_cut /= 0) THEN
        DO l=1,n3
           CALL shift_field (iim_file,jjm_file,i_cut,champ_file(:,:,l))
        ENDDO
     ENDIF
     ! interpolate onto the model grid
     DO l=1,n3
        DO j=1,jjm_g
           DO i=1,iim_g
              ij = i+(j-1)*iim_g
              champout_g(ij,l) = 0.
              DO m=1,ncorr(i,j)
                 champout_g(ij,l) = champout_g(ij,l) &
        &                        +champ_file(icorr(i,j,m),jcorr(i,j,m),l)
              ENDDO
              champout_g(ij,l) = champout_g(ij,l)/REAL(ncorr(i,j))
           ENDDO
        ENDDO
     ENDDO
!$     DO l=1,n3
!$        DO j=1,jjm_g
!$           WRITE(numout,*) j,(/ ( champout_g((j-1)*iim_g+i,l), i=1,iim_g ) /)
!$        ENDDO
!$     ENDDO
  ELSE
     ALLOCATE(champout_g(1,1))
  ENDIF
!!$  CALL scatter2D(champout_g,champout)
#ifndef CPP_PARA
  champout(:,:)=champout_g(:,:)
#else
  CALL scatter2D_rgen(champout_g,champout,n3)
#endif

!$  DO l=1,n3
!$     DO j=1,jjm
!$        WRITE(numout,*) j,(/ ( champout((j-1)*iim_g+i,l), i=1,iim_g ) /)
!$     ENDDO
!$  ENDDO
  !----------------------------
END SUBROUTINE weather_read
!-
!===
!-
!! ================================================================================================================================
!! SUBROUTINE   : weathgen_dump
!!
!>\BRIEF        This subroutine "dumps" data before writing an netcdf file.
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S):
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE weathgen_dump &
& (itau, dt_force, iim, jjm, nbindex, kindex, lrstwrite, &
&  swdown, rainf, snowf, tair, u, v, qair, pb, lwdown )
!---------------------------------------------------------------------
  IMPLICIT NONE
!-

  !! 0. Variables and parameters declaration

  !! 0.1 Input variables

  INTEGER,INTENT(IN)                     :: itau
  REAL,INTENT(IN)                        :: dt_force
  INTEGER,INTENT(IN)                     :: iim,jjm
  INTEGER,INTENT(IN)                     :: nbindex
  LOGICAL,INTENT(IN)                     :: lrstwrite
  INTEGER,DIMENSION(iim*jjm),INTENT(IN)  :: kindex
  REAL,DIMENSION(iim*jjm),INTENT(IN)     :: &
 &  swdown, rainf, snowf, tair, u, v, qair, pb, lwdown

  !! 0.4 Local variables

  INTEGER                 :: iret,ndim
  INTEGER,DIMENSION(4)    :: corner,edges
  REAL,DIMENSION(iim*jjm) :: var_gather

!_ ================================================================================================================================

!-
! time dimension
!-
  iret = NF90_PUT_VAR (dump_id,timestp_id,(/ REAL(itau) /), &
 &                     start=(/ itau /),count=(/ 1 /))
  iret = NF90_PUT_VAR (dump_id,time_id,(/ REAL(itau)*dt_force /), &
 &                     start=(/ itau /),count=(/ 1 /))
!-
! 2D variables: pas de dimension verticale
!-
  IF (gathered) THEN
    ndim = 2
    corner(1:2) = (/ 1, itau /)
    edges(1:2) = (/ nbindex, 1 /)
  ELSE
    ndim = 3
    corner(1:3) = (/ 1, 1, itau /)
    edges(1:3) = (/ iim, jjm, 1 /)
  ENDIF
!-
  CALL gather_weather (iim*jjm,nbindex,kindex,swdown, var_gather)
  iret = NF90_PUT_VAR (dump_id,soldownid, var_gather, &
 &         start=corner(1:ndim), count=edges(1:ndim))
  CALL gather_weather (iim*jjm,nbindex,kindex,rainf,  var_gather)
  iret = NF90_PUT_VAR (dump_id,rainfid,   var_gather, &
 &         start=corner(1:ndim), count=edges(1:ndim))
  CALL gather_weather (iim*jjm,nbindex,kindex,snowf,  var_gather)
  iret = NF90_PUT_VAR (dump_id,snowfid,   var_gather, &
 &         start=corner(1:ndim), count=edges(1:ndim))
  CALL gather_weather (iim*jjm,nbindex,kindex,pb,     var_gather)
  iret = NF90_PUT_VAR (dump_id,psolid,    var_gather, &
 &         start=corner(1:ndim), count=edges(1:ndim))
  CALL gather_weather (iim*jjm,nbindex,kindex,lwdown, var_gather)
  iret = NF90_PUT_VAR (dump_id,lwradid,   var_gather, &
 &         start=corner(1:ndim), count=edges(1:ndim))
!-
! 3D variables
!-
  IF (gathered) THEN
    ndim = 3
    corner(1:3) = (/ 1, 1, itau /)
    edges(1:3) = (/ nbindex, 1, 1 /)
  ELSE
    ndim = 4
    corner(1:4) = (/ 1, 1, 1, itau /)
    edges(1:4) = (/ iim, jjm, 1, 1 /)
  ENDIF
!-
  CALL gather_weather (iim*jjm,nbindex,kindex,u,    var_gather)
  iret = NF90_PUT_VAR (dump_id,uid,    var_gather, &
 &         start=corner(1:ndim), count=edges(1:ndim))
  CALL gather_weather (iim*jjm,nbindex,kindex,v,    var_gather)
  iret = NF90_PUT_VAR (dump_id,vid,    var_gather, &
 &         start=corner(1:ndim), count=edges(1:ndim))
  CALL gather_weather (iim*jjm,nbindex,kindex,tair, var_gather)
  iret = NF90_PUT_VAR (dump_id,tairid, var_gather, &
 &         start=corner(1:ndim), count=edges(1:ndim))
  CALL gather_weather (iim*jjm,nbindex,kindex,qair, var_gather)
  iret = NF90_PUT_VAR (dump_id,qairid, var_gather, &
 &         start=corner(1:ndim), count=edges(1:ndim))
!-
  IF (lrstwrite) THEN
    iret = NF90_CLOSE (dump_id)
  ENDIF
!---------------------------
END SUBROUTINE weathgen_dump
!-
!===
!-

!! ================================================================================================================================
!! SUBROUTINE   : gather_weather
!!
!>\BRIEF        This subroutine determines which points are not in the computational domain and
!! creates a mask for these points. 
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::var_out
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE gather_weather (iimjjm, nbindex, kindex, var_in, var_out)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-

  !! 0. Variables and parameters declaration

  !! 0.1 Input variables

  INTEGER,INTENT(IN)                     :: iimjjm,nbindex
  INTEGER,DIMENSION(iimjjm),INTENT(IN)   :: kindex
  REAL,DIMENSION(iimjjm),INTENT(IN)      :: var_in

  !! 0.2 Output variables  

  REAL,DIMENSION(iimjjm),INTENT(OUT)     :: var_out

  !! 0.4 Local variables

  INTEGER                                :: i
  LOGICAL,SAVE                           :: firstcall = .TRUE.
  INTEGER,SAVE                           :: nb_outside
  INTEGER,ALLOCATABLE,SAVE,DIMENSION(:)  :: outside

!_ ================================================================================================================================

  IF (firstcall) THEN
!---
!-- determine which points are not in the computational domain and
!-- create a mask for these points
!---
    firstcall = .FALSE.
  
    ALLOC_ERR=-1
    ALLOCATE(outside(iimjjm), STAT=ALLOC_ERR)
    IF (ALLOC_ERR/=0) THEN
      WRITE(numout,*) "ERROR IN ALLOCATION of outside : ",ALLOC_ERR
      STOP 
    ENDIF
    outside(:) = 0.
    nb_outside = 0
    DO i=1,iimjjm
      IF ( ALL( kindex(:) /= i ) ) THEN
        nb_outside = nb_outside+1
        outside(nb_outside) = i
      ENDIF
    ENDDO
  ENDIF
!-
  IF ( gathered ) THEN
    DO i=1,nbindex
      var_out(i) = var_in(kindex(i))
    ENDDO
  ELSE
    var_out(:) = var_in(:)
    DO i=1,nb_outside
      var_out(outside(i)) = undef
    ENDDO
  ENDIF
!--------------------
END SUBROUTINE gather_weather
!-
!===
!-

!! ================================================================================================================================
!! SUBROUTINE   : shift_shield
!!
!>\BRIEF        
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::champ
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE shift_field (im,jm,i_cut,champ)
!---------------------------------------------------------------------

  !! 0. Variables and parameters declaration

  !! 0.1 Input variables

  INTEGER,INTENT(IN)                  :: im,jm,i_cut

  !! 0.3 Modified variables

  REAL,DIMENSION(im,jm),INTENT(INOUT) :: champ

  !! 0.4 Local variables

  REAL,DIMENSION(im,jm)               :: champ_temp

!_ ================================================================================================================================

  IF ( (i_cut >= 1).AND.(i_cut <= im) ) THEN
    champ_temp(1:im-i_cut-1,:) = champ(i_cut:im,:)
    champ_temp(im-i_cut:im,:)  = champ(1:i_cut+1,:)
    champ(:,:) = champ_temp(:,:)
  ENDIF
!-------------------------
END SUBROUTINE shift_field
!-
!===
!-
!! ================================================================================================================================
!! SUBROUTINE   : weathgen_domain_size
!!
!>\BRIEF        This subroutine determines the size of the domain.
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::iim, ::jjm, ::limit_west, ::limit_east, ::limit_north, ::limit_south
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE weathgen_domain_size &
& (limit_west,limit_east,limit_north,limit_south, &
&  zonal_res,merid_res,iim,jjm)
!---------------------------------------------------------------------
  IMPLICIT NONE
!-
  !! 0. Variables and parameters declaration

  !! 0.1 Input variables

  REAL,INTENT(IN)     :: zonal_res,merid_res

  !! 0.2 Output variables 

  INTEGER,INTENT(OUT) :: iim,jjm

  !! 0.3 Modified variables

  REAL,INTENT(INOUT)  :: limit_west,limit_east,limit_north,limit_south

!_ ================================================================================================================================

  IF (limit_west > limit_east)  limit_east = limit_east+360.
!-
  IF (    (limit_west >= limit_east) &
 &    .OR.(limit_east > 360.) &
 &    .OR.(limit_west < -180.) &
 &    .OR.(limit_east-limit_west > 360.) ) THEN
    WRITE(numout,*) 'PROBLEME longitudes.'
    WRITE(numout,*) 'Limites Ouest, Est: ',limit_west,limit_east
    STOP
  ENDIF
!-
  IF (    (limit_south < -90.) &
 &    .OR.(limit_north > 90.) &
 &    .OR.(limit_south >= limit_north ) ) THEN
    WRITE(numout,*) 'PROBLEME latitudes.'
    WRITE(numout,*) 'Limites Nord, Sud: ',limit_north,limit_south
    STOP
  ENDIF
!-
  IF (    (zonal_res <= 0. ) &
 &    .OR.(zonal_res > limit_east-limit_west) ) THEN
    WRITE(numout,*) 'PROBLEME resolution zonale.'
    WRITE(numout,*) 'Limites Ouest, Est, Resolution: ', &
 &             limit_west,limit_east,zonal_res
    STOP
  ENDIF
!-
  IF (    (merid_res <= 0.) &
 &    .OR.(merid_res > limit_north-limit_south) ) THEN
    WRITE(numout,*) 'PROBLEME resolution meridionale.'
    WRITE(numout,*) 'Limites Nord, Sud, Resolution: ', &
 &             limit_north,limit_south,merid_res
    STOP
  ENDIF
!-
  iim = NINT(MAX((limit_east-limit_west)/zonal_res,1.))
  jjm = NINT(MAX((limit_north-limit_south)/merid_res,1.))
!-
  WRITE(numout,*) 'Domain size: iim, jjm = ', iim, jjm
!----------------------------------
END SUBROUTINE weathgen_domain_size
!-
!===
!-

!! ================================================================================================================================
!! FUNCTION     : tsat1 
!!
!>\BRIEF         This function computes the saturated temperature for vapor. 
!!
!! DESCRIPTION :  functions tsatl,tsati are used below so that Lowe's
!!                polyomial for liquid is used if t gt 273.16, or for ice if
!!                t lt 273.16. also impose range of validity for Lowe's polys.\n
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : tsat
!!
!! REFERENCE(S) : 
!! - Lowe PR, 1977, An Approximating Polynomial for the Computation of Saturation
!! Vapor Pressure, Journal of Applied Meteorology.16,101-103.
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

FUNCTION tsatl (t) RESULT (tsat)

  !! 0. Variables and parameters declaration
  
  !! 0.1 Input variables

  REAL,INTENT(IN) :: t

  !! 0.2 Result
  REAL            :: tsat

!_ ================================================================================================================================

  tsat = MIN(100.,MAX(t-zero_t,0.))
!-----------------
END FUNCTION tsatl
!-
!===
!-
!! ================================================================================================================================
!! FUNCTION     : tsati 
!!
!>\BRIEF         This function computes the saturated temperature for ice. 
!!
!! DESCRIPTION :  functions tsatl,tsati are used below so that Lowe's
!!                polyomial for liquid is used if t gt 273.16, or for ice if
!!                t lt 273.16. also impose range of validity for Lowe's polys.\n
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : tsat
!!
!! REFERENCE(S) : 
!! - Lowe PR, 1977, An Approximating Polynomial for the Computation of Saturation
!! Vapor Pressure, Journal of Applied Meteorology.16,101-103.
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

FUNCTION tsati (t) RESULT (tsat)

  !! 0. Variables and parameters declaration
  
  !! 0.1 Input variables

  REAL,INTENT(IN) :: t

  !! 0.2 Result
  REAL            :: tsat

!_ ================================================================================================================================

  tsat = MAX(-60.,MIN(t-zero_t,0.))
!-----------------
END FUNCTION tsati
!-
!===
!-

!! ================================================================================================================================
!! FUNCTION     : esat
!!
!>\BRIEF         This function computes specific humidity with the polynomials of Lowe.
!!
!! DESCRIPTION : The function esat is svp in n/m**2, with t in deg k.
!!               (100 * lowe's poly since 1 mb = 100 n/m**2.). \n
!!               Polynomials for svp(t), d(svp)/dt over water and ice 
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : esatout
!!
!! REFERENCE(S) : 
!! - Lowe PR, 1977, An Approximating Polynomial for the Computation of Saturation
!! Vapor Pressure, Journal of Applied Meteorology.16,101-103.
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

FUNCTION esat (t) RESULT (esatout)

  !! 0. Variables and parameters declaration

  REAL,PARAMETER :: &                                                    !! polynomials for svp(t), d(svp)/dt over water and ice are from
 & asat0 = 6.1078000,    asat1 = 4.4365185e-1, asat2 = 1.4289458e-2, &   !! lowe(1977),jam,16,101-103.
 & asat3 = 2.6506485e-4, asat4 = 3.0312404e-6, asat5 = 2.0340809e-8, &
 & asat6 = 6.1368209e-11, &
 & bsat0 = 6.1091780,    bsat1 = 5.0346990e-1, bsat2 = 1.8860134e-2, &
 & bsat3 = 4.1762237e-4, bsat4 = 5.8247203e-6, bsat5 = 4.8388032e-8, &
 & bsat6 = 1.8388269e-10

  !! 0.1 Input variables

  REAL,INTENT(IN) :: t

  !! 0.2 Result

  REAL             :: esatout

  !! 0.4 Local variables

  REAL             :: x

!_ ================================================================================================================================

  IF (t >= zero_t) THEN
    x = asat0
  ELSE
    x = bsat0
  ENDIF
!-
  esatout = 100.* &
    ( x &
     +tsatl(t)*(asat1+tsatl(t)*(asat2+tsatl(t)*(asat3 &
     +tsatl(t)*(asat4+tsatl(t)*(asat5+tsatl(t)* asat6))))) &
     +tsati(t)*(bsat1+tsati(t)*(bsat2+tsati(t)*(bsat3 &
     +tsati(t)*(bsat4+tsati(t)*(bsat5+tsati(t)* bsat6)))))  )
!----------------
END FUNCTION esat
!-
!===
!-
! ================================================================================================================================
!! FUNCTION     : qsat
!!
!>\BRIEF         This function computes the saturation specific humidity. 
!!
!! DESCRIPTION :  statement function qsat is saturation specific humidity,
!! with svp e and ambient pressure p in n/m**2. impose an upper
!! limit of 1 to avoid spurious values for very high svp
!! and/or small p
!!
!! RECENT CHANGE(S): None
!!
!! RETURN VALUE : qsatout
!!
!! REFERENCE(S) : 
!! - Lowe PR, 1977, An Approximating Polynomial for the Computation of Saturation
!! Vapor Pressure, Journal of Applied Meteorology.16,101-103.
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

FUNCTION qsat (e,p) RESULT (qsatout)

  !! 0. Variables and parameters declaration

  !! 0.1 Input variables

  REAL, INTENT(IN) :: e,p

  !! 0.2 Result
  REAL             :: qsatout

!_ ================================================================================================================================

  qsatout = 0.622*e/MAX(p-(1.0-0.622)*e,0.622*e)
!----------------
END FUNCTION qsat
!-
!===
!-

!! ================================================================================================================================
!! SUBROUTINE   : weathgen_qsat
!!
!>\BRIEF        This subroutine calculates the saturation vapor pressure with the polynomials of Lowe. 
!!
!! DESCRIPTION  : Vectorized version of functions esat and qsat.
!!                statement function esat is svp in n/m**2, with t in deg K.\n
!!               (100 * lowe's poly since 1 mb = 100 n/m**2.)\n
!!               Polynomials for svp(t), d(svp)/dt over water 
!!               and ice are from Lowe(1977),jam,16,101-103 \n
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::qsat
!!
!! REFERENCE(S) :
!! - Lowe PR, 1977, An Approximating Polynomial for the Computation of Saturation 
!! Vapor Pressure, Journal of Applied Meteorology.16,101-103.
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE weathgen_qsat (npoi,t,p,qsat)


  !! 0. Variables and parameters declaration

 !-
 ! polynomials for svp(t), d(svp)/dt over water and ice
 ! are from lowe(1977),ja
  REAL,PARAMETER :: &
 & asat0 = 6.1078000,    asat1 = 4.4365185e-1, asat2 = 1.4289458e-2, &   !! polynomials for svp(t), d(svp)/dt over water and ice
 & asat3 = 2.6506485e-4, asat4 = 3.0312404e-6, asat5 = 2.0340809e-8, &
 & asat6 = 6.1368209e-11, &
 & bsat0 = 6.1091780,    bsat1 = 5.0346990e-1, bsat2 = 1.8860134e-2, &
 & bsat3 = 4.1762237e-4, bsat4 = 5.8247203e-6, bsat5 = 4.8388032e-8, &
 & bsat6 = 1.8388269e-10

  !! 0.1 Input variables

  INTEGER,INTENT(IN)              :: npoi
  REAL,DIMENSION(npoi),INTENT(IN) :: t,p

  !! 0.2 Output variables 

  REAL,DIMENSION(npoi),INTENT(OUT):: qsat

  !! 0.4 Local variables

  REAL,DIMENSION(npoi) :: x, tl, ti, e

!_ ================================================================================================================================

  WHERE (t(:) > zero_t)
    x(:) = asat0
  ELSEWHERE
    x(:) = bsat0
  ENDWHERE
!-
  tl(:) = MIN(100.,MAX(t(:)-zero_t,0.))
  ti(:) = MAX(-60.,MIN(t(:)-zero_t,0.))
!-
  e(:) =  100.* &
    ( x(:) &
     +tl(:)*(asat1+tl(:)*(asat2+tl(:)*(asat3 &
     +tl(:)*(asat4+tl(:)*(asat5+tl(:)* asat6))))) &
     +ti(:)*(bsat1+ti(:)*(bsat2+ti(:)*(bsat3 &
     +ti(:)*(bsat4+ti(:)*(bsat5+ti(:)* bsat6)))))  )
!-
  qsat(:) = 0.622*e(:)/MAX(p(:)-(1.0-0.622)*e(:),0.622*e(:))
!---------------------------
END SUBROUTINE weathgen_qsat
!-
!===
!-

!! ================================================================================================================================
!! SUBROUTINE   : mask_c_o 
!!
!>\BRIEF        This subroutine computes a field which indicated if the point is terrestrial or oceanic.
!!
!! DESCRIPTION  : From a mask field, we make an indicator mask Earth/ocean. Earth : 1 ; ocean : 0.
!!
!! RECENT CHANGE(S): z.x.li,01/04/1994 : creation of the subroutine 
!!
!! MAIN OUTPUT VARIABLE(S): ::mask, ::ncorr, ::icorr
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE mask_c_o &
& (imdep, jmdep, xdata, ydata, mask_in, fcrit, &
&  imar, jmar, zonal_res, merid_res, n_agg, x, y, mask, &
&  ncorr, icorr, jcorr)

  !! 0. Variables and parameters declaration

  !! 0.1 Input variables

  INTEGER :: imdep,jmdep
  REAL :: xdata(imdep),ydata(jmdep)
  REAL :: mask_in(imdep,jmdep)
  REAL :: fcrit
  INTEGER :: imar,jmar
  REAL :: zonal_res,merid_res
  INTEGER :: n_agg
  REAL :: x(imar),y(jmar)

  !! 0.2 Output variables 

  REAL, INTENT(OUT) :: mask(imar,jmar)
  INTEGER :: ncorr(imar,jmar)
  INTEGER,DIMENSION(imar,jmar,n_agg) :: icorr,jcorr

  !! 0.4 Local variables 

  INTEGER i, j, ii, jj
  REAL a(imar),b(imar),c(jmar),d(jmar)
  INTEGER num_tot(imar,jmar), num_oce(imar,jmar)
  REAL,ALLOCATABLE :: distans(:)
  INTEGER ij_proche(1),i_proche,j_proche
!-
  INTEGER,DIMENSION(imar,jmar) :: ncorr_oce , ncorr_land
  INTEGER,DIMENSION(imar,jmar,n_agg) :: &
 &  icorr_oce, jcorr_oce , icorr_land, jcorr_land
!-
  INTEGER imdepp
  REAL,ALLOCATABLE :: xdatap(:)
  REAL,ALLOCATABLE :: mask_inp(:,:)
  LOGICAL :: extend

!_ ================================================================================================================================

  ncorr(:,:) = 0
  icorr(:,:,:) = -1; jcorr(:,:,:) = -1
  ncorr_land(:,:) = 0
  icorr_land(:,:,:) = -1; jcorr_land(:,:,:) = -1
  ncorr_oce(:,:) = 0
  icorr_oce(:,:,:) = -1; jcorr_oce(:,:,:) = -1
! do we have to extend the domain (-x...-x+360)?
  IF ( xdata(1)+360.-xdata(imdep) > 0.01 ) THEN
    extend = .TRUE.
    imdepp = imdep+1
  ELSE
    extend = .FALSE.
    imdepp = imdep
  ENDIF
!-

  ALLOC_ERR=-1
  ALLOCATE(xdatap(imdepp), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of xdatap : ",ALLOC_ERR
    STOP 
  ENDIF

  ALLOC_ERR=-1
  ALLOCATE(mask_inp(imdepp,jmdep), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of mask_inp : ",ALLOC_ERR
    STOP 
  ENDIF
!-
  xdatap(1:imdep) = xdata(1:imdep)
  mask_inp(1:imdep,:) = mask_in(1:imdep,:)
!-
  IF (extend) THEN
    xdatap(imdepp) = xdatap(1)+360.
    mask_inp(imdepp,:) = mask_inp(1,:)
  ENDIF
!-

  ALLOC_ERR=-1
  ALLOCATE(distans(imdepp*jmdep), STAT=ALLOC_ERR)
  IF (ALLOC_ERR/=0) THEN
    WRITE(numout,*) "ERROR IN ALLOCATION of distans : ",ALLOC_ERR
    STOP 
  ENDIF
! Definition des limites des boites de la grille d'arrivee.
  IF (imar > 1) THEN
    a(1) = x(1)-(x(2)-x(1))/2.0
    b(1) = (x(1)+x(2))/2.0
    DO i=2,imar-1
      a(i) = b(i-1)
      b(i) = (x(i)+x(i+1))/2.0
    ENDDO
    a(imar) = b(imar-1)
    b(imar) = x(imar)+(x(imar)-x(imar-1))/2.0
  ELSE
    a(1) = x(1)-zonal_res/2.
    b(1) = x(1)+zonal_res/2.
  ENDIF
!-
  IF (jmar > 1) THEN
    c(1) = y(1)-(y(2)-y(1))/2.0
    d(1) = (y(1)+y(2))/2.0
    DO j=2,jmar-1
      c(j) = d(j-1)
      d(j) = (y(j)+y(j+1))/2.0
    ENDDO
    c(jmar) = d(jmar-1)
    d(jmar) = y(jmar)+(y(jmar)-y(jmar-1))/2.0
  ELSE
    c(1) = y(1)-merid_res/2.
    d(1) = y(1)+merid_res/2.
  ENDIF
!-
  num_oce(1:imar,1:jmar) = 0
  num_tot(1:imar,1:jmar) = 0
!-
!  .....  Modif  P. Le Van ( 23/08/95 )  ....
!-
  DO ii=1,imar
    DO jj=1,jmar
      DO i=1,imdepp
        IF (    (     (xdatap(i)-a(ii) >= 1.e-5) &
 &               .AND.(xdatap(i)-b(ii) <= 1.e-5) ) &
 &          .OR.(     (xdatap(i)-a(ii) <= 1.e-5) &
 &               .AND.(xdatap(i)-b(ii) >= 1.e-5) ) ) THEN
          DO j=1,jmdep
            IF (    (     (ydata(j)-c(jj) >= 1.e-5) &
 &                   .AND.(ydata(j)-d(jj) <= 1.e-5) ) &
 &              .OR.(     (ydata(j)-c(jj) <= 1.e-5) &
 &                   .AND.(ydata(j)-d(jj) >= 1.e-5) ) ) THEN
              num_tot(ii,jj) = num_tot(ii,jj)+1
              IF (mask_inp(i,j) < 0.5) THEN
                num_oce(ii,jj) = num_oce(ii,jj)+1
!-------------- on a trouve un point oceanique. On le memorise.
                ncorr_oce(ii,jj) = ncorr_oce(ii,jj)+1
                IF ((i == imdepp).AND.extend) THEN
                  icorr_oce(ii,jj,ncorr_oce(ii,jj)) = 1
                ELSE
                  icorr_oce(ii,jj,ncorr_oce(ii,jj)) = i
                ENDIF
                jcorr_oce(ii,jj,ncorr_oce(ii,jj)) = j
              ELSE
!-------------- on a trouve un point continental. On le memorise.
                ncorr_land(ii,jj) = ncorr_land(ii,jj)+1
                IF ((i == imdepp).AND.extend) THEN
                  icorr_land(ii,jj,ncorr_land(ii,jj)) = 1
                ELSE
                  icorr_land(ii,jj,ncorr_land(ii,jj)) = i
                ENDIF
                jcorr_land(ii,jj,ncorr_land(ii,jj)) = j
              ENDIF
            ENDIF
          ENDDO
        ENDIF
      ENDDO
    ENDDO
  ENDDO
!-
  DO i=1,imar
    DO j=1,jmar
      IF (num_tot(i,j) > 0) THEN
        IF (    (     (num_oce(i,j) == 0) &
 &               .AND.(num_tot(i,j) > 0) ) &
 &          .OR.(     (num_oce(i,j) > 0) &
 &               .AND.(   REAL(num_oce(i,j)) &
 &                     <= REAL(num_tot(i,j))*(1.-fcrit) ) ) ) THEN
          mask(i,j) = 1.
          ncorr(i,j) = ncorr_land(i,j)
          icorr(i,j,:) = icorr_land(i,j,:)
          jcorr(i,j,:) = jcorr_land(i,j,:)
        ELSE
          mask(i,j) = 0.
          ncorr(i,j) = ncorr_oce(i,j)
          icorr(i,j,:) = icorr_oce(i,j,:)
          jcorr(i,j,:) = jcorr_oce(i,j,:)
        ENDIF
      ELSE
        CALL dist_sphe(x(i),y(j),xdatap,ydata,imdepp,jmdep,distans)
        ij_proche(:) = MINLOC(distans)
        j_proche = (ij_proche(1)-1)/imdepp+1
        i_proche = ij_proche(1)-(j_proche-1)*imdepp
        mask(i,j) = mask_inp(i_proche,j_proche)
        IF ( (i_proche == imdepp).AND.extend)  i_proche=1
        ncorr(i,j) = 1
        icorr(i,j,1) = i_proche
        jcorr(i,j,1) = j_proche
      ENDIF
    ENDDO
  ENDDO
!----------------------
END SUBROUTINE mask_c_o
!-
!===
!-
!! ================================================================================================================================
!! SUBROUTINE   : dist_sphe
!!
!>\BRIEF        This subroutine computes the minimum distance between two points on Earth
!! according the great circle. \n
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): Laurent Li, 30/12/1996 : creation of this subroutine
!!
!! MAIN OUTPUT VARIABLE(S): ::distance
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE dist_sphe (rf_lon,rf_lat,rlon,rlat,im,jm,distance)


  !! 0. Variables and parameters declaration

  REAL radius                !! Earth radius (m)
  PARAMETER (radius=6371229.)

  !! 0.1 Input variables

  INTEGER im, jm    !! dimensions
  REAL rf_lon       !! longitude of the referenced point (degrees)
  REAL rf_lat       !! latitude of the referenced point (degrees)
  REAL rlon(im), rlat(jm) !! longitude and latitude of points

  !! 0.2 Output variables

  REAL distance(im,jm) !! distances in meters (m)

  !! 0.4 Local variables

  REAL rlon1, rlat1
  REAL rlon2, rlat2
  REAL dist
  REAL pa, pb, p, pi
  INTEGER i,j

!_ ================================================================================================================================

  pi = 4.0*ATAN(1.0)
!-
  DO j=1,jm
    DO i=1,im
      rlon1=rf_lon
      rlat1=rf_lat
      rlon2=rlon(i)
      rlat2=rlat(j)
      pa = pi/2.0-rlat1*pi/180.0 ! dist. entre pole n et point a
      pb = pi/2.0-rlat2*pi/180.0 ! dist. entre pole n et point b
!-----
      p = (rlon1-rlon2)*pi/180.0 ! angle entre a et b (leurs meridiens)
!-----
      dist = ACOS( COS(pa)*COS(pb)+SIN(pa)*SIN(pb)*COS(p))
      dist = radius*dist
      distance(i,j) = dist
    ENDDO
  ENDDO
!-----------------------
END SUBROUTINE dist_sphe
!-
!===
!-

!! ================================================================================================================================
!! SUBROUTINE   : permute
!!
!>\BRIEF         This subroutine initializes an array of size n by random integers between 1 and n.  
!!
!! DESCRIPTION  : None
!!
!! RECENT CHANGE(S): None
!!
!! MAIN OUTPUT VARIABLE(S): ::ordre
!!
!! REFERENCE(S) :
!!
!! FLOWCHART    : None
!! \n
!_ ================================================================================================================================

SUBROUTINE permute (n,ordre)
!---------------------------------------------------------------------

  !! 0. Variables and parameters declaration 

  !! 0.1 Input variables

  INTEGER,INTENT(IN) :: n   !! size of the array

  !! 0.2 Output variables 

  INTEGER,DIMENSION(n),INTENT(OUT) :: ordre

  !! 0.4 Local variables

  INTEGER,DIMENSION(n) :: restant
  INTEGER :: ipique, i, n_rest
  REAL    :: rndnum

!_ ================================================================================================================================

  n_rest = n
  restant(:) = (/ (i, i=1,n) /)
!-
  DO i=1,n
    CALL random_number (rndnum)
    ipique = INT(rndnum*n_rest)+1
    ordre(i) = restant(ipique)
    restant(ipique:n_rest-1) = restant(ipique+1:n_rest)
    n_rest = n_rest-1
  ENDDO
!---------------------
END SUBROUTINE permute
!-
!===
!-----------------
END MODULE weather