trunk/dyn3d/guide.f

MODULE guide_m

  ! From dyn3d/guide.F, version 1.3 2005/05/25 13:10:09
  ! and dyn3d/guide.h, version 1.1.1.1 2004/05/19 12:53:06

  IMPLICIT NONE

  REAL tau_min_u, tau_max_u
  REAL tau_min_v, tau_max_v
  REAL tau_min_t, tau_max_t
  REAL tau_min_q, tau_max_q
  REAL tau_min_p, tau_max_p
  REAL aire_min, aire_max


  LOGICAL guide_u, guide_v, guide_t, guide_q, guide_p
  LOGICAL ncep, ini_anal
  INTEGER online

CONTAINS

  SUBROUTINE guide(itau, ucov, vcov, teta, q, masse, ps)

    ! Author: F.Hourdin

    USE dimens_m, ONLY : jjm, llm
    USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1, llmp1
    USE comconst, ONLY : cpp, daysec, dtvr, kappa, pi
    USE comvert, ONLY : ap, bp, preff, presnivs
    USE conf_gcm_m, ONLY : day_step, iperiod
    USE comgeom, ONLY : aire, rlatu, rlonv
    USE serre, ONLY : clat, clon
    USE q_sat_m, ONLY : q_sat
    USE exner_hyb_m, ONLY : exner_hyb
    USE inigrads_m, ONLY : inigrads
    use netcdf, only: nf90_nowrite, nf90_open, nf90_close
    use tau2alpha_m, only: tau2alpha

    INCLUDE 'netcdf.inc'

    !   variables dynamiques
    REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
    REAL, intent(inout):: teta(ip1jmp1, llm) ! temperature potentielle 
    REAL q(ip1jmp1, llm) ! temperature potentielle 
    REAL ps(ip1jmp1) ! pression  au sol
    REAL masse(ip1jmp1, llm) ! masse d'air

    !   common passe pour des sorties
    REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
    COMMON /comdxdy/dxdys, dxdyu, dxdyv

    !   variables dynamiques pour les reanalyses.
    REAL ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas
    REAL tetarea1(ip1jmp1, llm) ! temp pot  reales
    REAL qrea1(ip1jmp1, llm) ! temp pot  reales
    REAL psrea1(ip1jmp1) ! ps
    REAL ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas
    REAL tetarea2(ip1jmp1, llm) ! temp pot  reales
    REAL qrea2(ip1jmp1, llm) ! temp pot  reales
    REAL masserea2(ip1jmp1, llm) ! masse
    REAL psrea2(ip1jmp1) ! ps

    REAL alpha_q(ip1jmp1)
    REAL alpha_t(ip1jmp1), alpha_p(ip1jmp1)
    REAL alpha_u(ip1jmp1), alpha_v(ip1jm)
    REAL dday_step, toto, reste, itau_test
    INTEGER step_rea, count_no_rea

    INTEGER ilon, ilat
    REAL factt, ztau(ip1jmp1)

    INTEGER, INTENT (IN) :: itau
    INTEGER ij, l
    INTEGER ncidpl, varidpl, nlev, status
    INTEGER rcod, rid
    REAL ditau, tau, a
    SAVE nlev

    !  TEST SUR QSAT
    REAL p(ip1jmp1, llmp1), pk(ip1jmp1, llm), pks(ip1jmp1)
    REAL pkf(ip1jmp1, llm)
    REAL pres(ip1jmp1, llm)

    REAL qsat(ip1jmp1, llm)
    REAL unskap
    REAL tnat(ip1jmp1, llm)

    LOGICAL:: first = .TRUE.

    SAVE ucovrea1, vcovrea1, tetarea1, psrea1, qrea1
    SAVE ucovrea2, vcovrea2, tetarea2, masserea2, psrea2, qrea2

    SAVE alpha_t, alpha_q, alpha_u, alpha_v, alpha_p, itau_test
    SAVE step_rea, count_no_rea

    CHARACTER (10) file
    INTEGER igrads
    REAL dtgrads
    SAVE igrads, dtgrads
    DATA igrads, dtgrads/2, 100./

    !-----------------------------------------------------------------------

    PRINT *, 'Call sequence information: guide'

    ! calcul de l'humidite saturante

    forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
    CALL massdair(p, masse)
    PRINT *, 'OK1'
    CALL exner_hyb(ps, p, pks, pk, pkf)
    PRINT *, 'OK2'
    tnat(:, :) = pk(:, :)*teta(:, :)/cpp
    PRINT *, 'OK3'
    unskap = 1./kappa
    pres(:, :) = preff*(pk(:, :)/cpp)**unskap
    PRINT *, 'OK4'
    qsat = q_sat(tnat, pres)

    !   initialisations pour la lecture des reanalyses.
    !    alpha determine la part des injections de donnees a chaque etape
    !    alpha=1 signifie pas d'injection
    !    alpha=0 signifie injection totale

    PRINT *, 'ONLINE=', online
    IF (online==-1) THEN
       RETURN
    END IF

    IF (first) THEN

       PRINT *, 'initialisation du guide '
       CALL conf_guide
       PRINT *, 'apres conf_guide'

       file = 'guide'
       CALL inigrads(igrads, rlonv, 180./pi, -180., 180., rlatu, -90., 90., &
            180./pi, presnivs, 1., dtgrads, file, 'dyn_zon ')

       PRINT *, &
            '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'

       IF (online==-1) RETURN
       IF (online==1) THEN

          !  Constantes de temps de rappel en jour
          !  0.1 c'est en gros 2h30. 
          !  1e10  est une constante infinie donc en gros pas de guidage

          !   coordonnees du centre du zoom
          CALL coordij(clon, clat, ilon, ilat)
          !   aire de la maille au centre du zoom
          aire_min = aire(ilon+(ilat-1)*iip1)
          !   aire maximale de la maille
          aire_max = 0.
          DO ij = 1, ip1jmp1
             aire_max = max(aire_max, aire(ij))
          END DO
          !  factt = pas de temps en fraction de jour
          factt = dtvr*iperiod/daysec

          CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
          CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
          CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)
          CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)
          CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)

          CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
          CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U   ')
          CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T   ')

          !   Cas ou on force exactement par les variables analysees
       ELSE
          alpha_t = 0.
          alpha_u = 0.
          alpha_v = 0.
          alpha_p = 0.
          !           physic=.false.
       END IF

       itau_test = 1001
       step_rea = 1
       count_no_rea = 0
       ncidpl = -99

       !    itau_test    montre si l'importation a deja ete faite au rang itau
       ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
       if (guide_u) then
          if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
       endif

       if (guide_v) then
          if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
       endif

       if (guide_T) then
          if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
       endif

       if (guide_Q) then
          if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
       endif

       IF (ncep) THEN
          status = nf_inq_dimid(ncidpl, 'LEVEL', rid)
       ELSE
          status = nf_inq_dimid(ncidpl, 'PRESSURE', rid)
       END IF
       status = nf_inq_dimlen(ncidpl, rid, nlev)
       PRINT *, 'nlev', nlev
       rcod = nf90_close(ncidpl)
       !   Lecture du premier etat des reanalyses.
       CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
            masserea2, psrea2, 1, nlev)
       qrea2(:, :) = max(qrea2(:, :), 0.1)


       !   Debut de l'integration temporelle:
    END IF ! first

    ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:

    ditau = real(itau)
    dday_step = real(day_step)
    WRITE (*, *) 'ditau, dday_step'
    WRITE (*, *) ditau, dday_step
    toto = 4*ditau/dday_step
    reste = toto - aint(toto)

    IF (reste==0.) THEN
       IF (itau_test==itau) THEN
          WRITE (*, *) 'deuxieme passage de advreel a itau=', itau
          STOP
       ELSE
          vcovrea1(:, :) = vcovrea2(:, :)
          ucovrea1(:, :) = ucovrea2(:, :)
          tetarea1(:, :) = tetarea2(:, :)
          qrea1(:, :) = qrea2(:, :)

          PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
               count_no_rea, ' non lectures'
          step_rea = step_rea + 1
          itau_test = itau
          CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &
               qrea2, masserea2, psrea2, 1, nlev)
          qrea2(:, :) = max(qrea2(:, :), 0.1)
          factt = dtvr*iperiod/daysec
          ztau(:) = factt/max(alpha_t(:), 1.E-10)
          CALL wrgrads(igrads, 1, aire, 'aire      ', 'aire      ')
          CALL wrgrads(igrads, 1, dxdys, 'dxdy      ', 'dxdy      ')
          CALL wrgrads(igrads, 1, alpha_u, 'au        ', 'au        ')
          CALL wrgrads(igrads, 1, alpha_t, 'at        ', 'at        ')
          CALL wrgrads(igrads, 1, ztau, 'taut      ', 'taut      ')
          CALL wrgrads(igrads, llm, ucov, 'u         ', 'u         ')
          CALL wrgrads(igrads, llm, ucovrea2, 'ua        ', 'ua        ')
          CALL wrgrads(igrads, llm, teta, 'T         ', 'T         ')
          CALL wrgrads(igrads, llm, tetarea2, 'Ta        ', 'Ta        ')
          CALL wrgrads(igrads, llm, qrea2, 'Qa        ', 'Qa        ')
          CALL wrgrads(igrads, llm, q, 'Q         ', 'Q         ')

          CALL wrgrads(igrads, llm, qsat, 'QSAT      ', 'QSAT      ')

       END IF
    ELSE
       count_no_rea = count_no_rea + 1
    END IF

    !   Guidage
    !    x_gcm = a * x_gcm + (1-a) * x_reanalyses

    IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'

    ditau = real(itau)
    dday_step = real(day_step)


    tau = 4*ditau/dday_step
    tau = tau - aint(tau)

    !  ucov
    IF (guide_u) THEN
       DO l = 1, llm
          DO ij = 1, ip1jmp1
             a = (1.-tau)*ucovrea1(ij, l) + tau*ucovrea2(ij, l)
             ucov(ij, l) = (1.-alpha_u(ij))*ucov(ij, l) + alpha_u(ij)*a
             IF (first .AND. ini_anal) ucov(ij, l) = a
          END DO
       END DO
    END IF

    IF (guide_t) THEN
       DO l = 1, llm
          DO ij = 1, ip1jmp1
             a = (1.-tau)*tetarea1(ij, l) + tau*tetarea2(ij, l)
             teta(ij, l) = (1.-alpha_t(ij))*teta(ij, l) + alpha_t(ij)*a
             IF (first .AND. ini_anal) teta(ij, l) = a
          END DO
       END DO
    END IF

    !  P
    IF (guide_p) THEN
       DO ij = 1, ip1jmp1
          a = (1.-tau)*psrea1(ij) + tau*psrea2(ij)
          ps(ij) = (1.-alpha_p(ij))*ps(ij) + alpha_p(ij)*a
          IF (first .AND. ini_anal) ps(ij) = a
       END DO
       forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
       CALL massdair(p, masse)
    END IF


    !  q
    IF (guide_q) THEN
       DO l = 1, llm
          DO ij = 1, ip1jmp1
             a = (1.-tau)*qrea1(ij, l) + tau*qrea2(ij, l)
             !   hum relative en % -> hum specif
             a = qsat(ij, l)*a*0.01
             q(ij, l) = (1.-alpha_q(ij))*q(ij, l) + alpha_q(ij)*a
             IF (first .AND. ini_anal) q(ij, l) = a
          END DO
       END DO
    END IF

    ! vcov
    IF (guide_v) THEN
       DO l = 1, llm
          DO ij = 1, ip1jm
             a = (1.-tau)*vcovrea1(ij, l) + tau*vcovrea2(ij, l)
             vcov(ij, l) = (1.-alpha_v(ij))*vcov(ij, l) + alpha_v(ij)*a
             IF (first .AND. ini_anal) vcov(ij, l) = a
          END DO
          IF (first .AND. ini_anal) vcov(ij, l) = a
       END DO
    END IF

    first = .FALSE.

  END SUBROUTINE guide

END MODULE guide_m
1	guez	20	MODULE guide_m
2	guez	3
3	guez	29	! From dyn3d/guide.F, version 1.3 2005/05/25 13:10:09
4			! and dyn3d/guide.h, version 1.1.1.1 2004/05/19 12:53:06
5	guez	3
6	guez	37	IMPLICIT NONE
7
8	guez	36	REAL tau_min_u, tau_max_u
9			REAL tau_min_v, tau_max_v
10			REAL tau_min_t, tau_max_t
11			REAL tau_min_q, tau_max_q
12			REAL tau_min_p, tau_max_p
13			REAL aire_min, aire_max
14	guez	3
15
16	guez	36	LOGICAL guide_u, guide_v, guide_t, guide_q, guide_p
17			LOGICAL ncep, ini_anal
18			INTEGER online
19	guez	3
20	guez	20	CONTAINS
21	guez	3
22	guez	29	SUBROUTINE guide(itau, ucov, vcov, teta, q, masse, ps)
23	guez	3
24	guez	29	! Author: F.Hourdin
25	guez	3
26	guez	29	USE dimens_m, ONLY : jjm, llm
27			USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1, llmp1
28			USE comconst, ONLY : cpp, daysec, dtvr, kappa, pi
29			USE comvert, ONLY : ap, bp, preff, presnivs
30			USE conf_gcm_m, ONLY : day_step, iperiod
31			USE comgeom, ONLY : aire, rlatu, rlonv
32			USE serre, ONLY : clat, clon
33			USE q_sat_m, ONLY : q_sat
34			USE exner_hyb_m, ONLY : exner_hyb
35			USE inigrads_m, ONLY : inigrads
36			use netcdf, only: nf90_nowrite, nf90_open, nf90_close
37	guez	37	use tau2alpha_m, only: tau2alpha
38	guez	3
39	guez	29	INCLUDE 'netcdf.inc'
40	guez	3
41	guez	29	! variables dynamiques
42	guez	36	REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
43	guez	29	REAL, intent(inout):: teta(ip1jmp1, llm) ! temperature potentielle
44	guez	36	REAL q(ip1jmp1, llm) ! temperature potentielle
45			REAL ps(ip1jmp1) ! pression au sol
46			REAL masse(ip1jmp1, llm) ! masse d'air
47	guez	3
48	guez	29	! common passe pour des sorties
49	guez	36	REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
50	guez	29	COMMON /comdxdy/dxdys, dxdyu, dxdyv
51	guez	3
52	guez	29	! variables dynamiques pour les reanalyses.
53	guez	36	REAL ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas
54			REAL tetarea1(ip1jmp1, llm) ! temp pot reales
55			REAL qrea1(ip1jmp1, llm) ! temp pot reales
56			REAL psrea1(ip1jmp1) ! ps
57			REAL ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas
58			REAL tetarea2(ip1jmp1, llm) ! temp pot reales
59			REAL qrea2(ip1jmp1, llm) ! temp pot reales
60			REAL masserea2(ip1jmp1, llm) ! masse
61			REAL psrea2(ip1jmp1) ! ps
62	guez	3
63	guez	36	REAL alpha_q(ip1jmp1)
64			REAL alpha_t(ip1jmp1), alpha_p(ip1jmp1)
65			REAL alpha_u(ip1jmp1), alpha_v(ip1jm)
66			REAL dday_step, toto, reste, itau_test
67			INTEGER step_rea, count_no_rea
68	guez	3
69	guez	36	INTEGER ilon, ilat
70			REAL factt, ztau(ip1jmp1)
71	guez	3
72	guez	29	INTEGER, INTENT (IN) :: itau
73	guez	36	INTEGER ij, l
74			INTEGER ncidpl, varidpl, nlev, status
75			INTEGER rcod, rid
76			REAL ditau, tau, a
77	guez	29	SAVE nlev
78	guez	3
79	guez	29	! TEST SUR QSAT
80	guez	36	REAL p(ip1jmp1, llmp1), pk(ip1jmp1, llm), pks(ip1jmp1)
81			REAL pkf(ip1jmp1, llm)
82			REAL pres(ip1jmp1, llm)
83	guez	3
84	guez	36	REAL qsat(ip1jmp1, llm)
85			REAL unskap
86			REAL tnat(ip1jmp1, llm)
87	guez	3
88	guez	37	LOGICAL:: first = .TRUE.
89	guez	3
90	guez	29	SAVE ucovrea1, vcovrea1, tetarea1, psrea1, qrea1
91			SAVE ucovrea2, vcovrea2, tetarea2, masserea2, psrea2, qrea2
92	guez	3
93	guez	29	SAVE alpha_t, alpha_q, alpha_u, alpha_v, alpha_p, itau_test
94			SAVE step_rea, count_no_rea
95	guez	3
96	guez	36	CHARACTER (10) file
97			INTEGER igrads
98			REAL dtgrads
99	guez	29	SAVE igrads, dtgrads
100			DATA igrads, dtgrads/2, 100./
101	guez	3
102	guez	29	!-----------------------------------------------------------------------
103	guez	3
104	guez	29	PRINT *, 'Call sequence information: guide'
105	guez	3
106	guez	29	! calcul de l'humidite saturante
107	guez	3
108	guez	37	forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
109	guez	29	CALL massdair(p, masse)
110			PRINT *, 'OK1'
111			CALL exner_hyb(ps, p, pks, pk, pkf)
112			PRINT *, 'OK2'
113			tnat(:, :) = pk(:, :)*teta(:, :)/cpp
114			PRINT *, 'OK3'
115			unskap = 1./kappa
116			pres(:, :) = preff(pk(:, :)/cpp)*unskap
117			PRINT *, 'OK4'
118			qsat = q_sat(tnat, pres)
119	guez	3
120	guez	29	! initialisations pour la lecture des reanalyses.
121			! alpha determine la part des injections de donnees a chaque etape
122			! alpha=1 signifie pas d'injection
123			! alpha=0 signifie injection totale
124	guez	3
125	guez	29	PRINT *, 'ONLINE=', online
126			IF (online==-1) THEN
127			RETURN
128			END IF
129	guez	3
130	guez	29	IF (first) THEN
131	guez	3
132	guez	29	PRINT *, 'initialisation du guide '
133			CALL conf_guide
134			PRINT *, 'apres conf_guide'
135	guez	3
136	guez	29	file = 'guide'
137			CALL inigrads(igrads, rlonv, 180./pi, -180., 180., rlatu, -90., 90., &
138			180./pi, presnivs, 1., dtgrads, file, 'dyn_zon ')
139	guez	3
140	guez	29	PRINT *, &
141			'1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
142	guez	3
143	guez	29	IF (online==-1) RETURN
144			IF (online==1) THEN
145	guez	3
146	guez	29	! Constantes de temps de rappel en jour
147			! 0.1 c'est en gros 2h30.
148			! 1e10 est une constante infinie donc en gros pas de guidage
149	guez	3
150	guez	29	! coordonnees du centre du zoom
151			CALL coordij(clon, clat, ilon, ilat)
152			! aire de la maille au centre du zoom
153			aire_min = aire(ilon+(ilat-1)*iip1)
154			! aire maximale de la maille
155			aire_max = 0.
156			DO ij = 1, ip1jmp1
157			aire_max = max(aire_max, aire(ij))
158			END DO
159			! factt = pas de temps en fraction de jour
160			factt = dtvr*iperiod/daysec
161	guez	3
162	guez	29	CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
163			CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
164			CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)
165			CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)
166			CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)
167	guez	3
168	guez	29	CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
169			CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')
170			CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')
171	guez	3
172	guez	29	! Cas ou on force exactement par les variables analysees
173			ELSE
174			alpha_t = 0.
175			alpha_u = 0.
176			alpha_v = 0.
177			alpha_p = 0.
178			! physic=.false.
179			END IF
180	guez	3
181	guez	29	itau_test = 1001
182			step_rea = 1
183			count_no_rea = 0
184			ncidpl = -99
185	guez	3
186	guez	29	! itau_test montre si l'importation a deja ete faite au rang itau
187			! lecture d'un fichier netcdf pour determiner le nombre de niveaux
188			if (guide_u) then
189			if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
190			endif
191	guez	3
192	guez	29	if (guide_v) then
193			if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
194			endif
195	guez	3
196	guez	29	if (guide_T) then
197			if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
198			endif
199	guez	3
200	guez	29	if (guide_Q) then
201			if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
202			endif
203	guez	3
204	guez	29	IF (ncep) THEN
205			status = nf_inq_dimid(ncidpl, 'LEVEL', rid)
206			ELSE
207			status = nf_inq_dimid(ncidpl, 'PRESSURE', rid)
208			END IF
209			status = nf_inq_dimlen(ncidpl, rid, nlev)
210			PRINT *, 'nlev', nlev
211			rcod = nf90_close(ncidpl)
212			! Lecture du premier etat des reanalyses.
213			CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
214			masserea2, psrea2, 1, nlev)
215			qrea2(:, :) = max(qrea2(:, :), 0.1)
216	guez	3
217
218	guez	29	! Debut de l'integration temporelle:
219			END IF ! first
220	guez	3
221	guez	29	! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
222	guez	3
223	guez	29	ditau = real(itau)
224			dday_step = real(day_step)
225			WRITE (, ) 'ditau, dday_step'
226			WRITE (, ) ditau, dday_step
227			toto = 4*ditau/dday_step
228			reste = toto - aint(toto)
229	guez	3
230	guez	29	IF (reste==0.) THEN
231			IF (itau_test==itau) THEN
232			WRITE (, ) 'deuxieme passage de advreel a itau=', itau
233			STOP
234			ELSE
235			vcovrea1(:, :) = vcovrea2(:, :)
236			ucovrea1(:, :) = ucovrea2(:, :)
237			tetarea1(:, :) = tetarea2(:, :)
238			qrea1(:, :) = qrea2(:, :)
239	guez	3
240	guez	29	PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
241			count_no_rea, ' non lectures'
242			step_rea = step_rea + 1
243			itau_test = itau
244			CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &
245			qrea2, masserea2, psrea2, 1, nlev)
246			qrea2(:, :) = max(qrea2(:, :), 0.1)
247			factt = dtvr*iperiod/daysec
248			ztau(:) = factt/max(alpha_t(:), 1.E-10)
249			CALL wrgrads(igrads, 1, aire, 'aire ', 'aire ')
250			CALL wrgrads(igrads, 1, dxdys, 'dxdy ', 'dxdy ')
251			CALL wrgrads(igrads, 1, alpha_u, 'au ', 'au ')
252			CALL wrgrads(igrads, 1, alpha_t, 'at ', 'at ')
253			CALL wrgrads(igrads, 1, ztau, 'taut ', 'taut ')
254			CALL wrgrads(igrads, llm, ucov, 'u ', 'u ')
255			CALL wrgrads(igrads, llm, ucovrea2, 'ua ', 'ua ')
256			CALL wrgrads(igrads, llm, teta, 'T ', 'T ')
257			CALL wrgrads(igrads, llm, tetarea2, 'Ta ', 'Ta ')
258			CALL wrgrads(igrads, llm, qrea2, 'Qa ', 'Qa ')
259			CALL wrgrads(igrads, llm, q, 'Q ', 'Q ')
260	guez	3
261	guez	29	CALL wrgrads(igrads, llm, qsat, 'QSAT ', 'QSAT ')
262	guez	3
263	guez	29	END IF
264			ELSE
265			count_no_rea = count_no_rea + 1
266			END IF
267	guez	3
268	guez	29	! Guidage
269			! x_gcm = a * x_gcm + (1-a) * x_reanalyses
270	guez	3
271	guez	29	IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'
272	guez	3
273	guez	29	ditau = real(itau)
274			dday_step = real(day_step)
275	guez	3
276
277	guez	29	tau = 4*ditau/dday_step
278			tau = tau - aint(tau)
279	guez	3
280	guez	29	! ucov
281			IF (guide_u) THEN
282			DO l = 1, llm
283			DO ij = 1, ip1jmp1
284			a = (1.-tau)ucovrea1(ij, l) + tauucovrea2(ij, l)
285			ucov(ij, l) = (1.-alpha_u(ij))ucov(ij, l) + alpha_u(ij)a
286			IF (first .AND. ini_anal) ucov(ij, l) = a
287			END DO
288			END DO
289			END IF
290	guez	3
291	guez	29	IF (guide_t) THEN
292			DO l = 1, llm
293			DO ij = 1, ip1jmp1
294			a = (1.-tau)tetarea1(ij, l) + tautetarea2(ij, l)
295			teta(ij, l) = (1.-alpha_t(ij))teta(ij, l) + alpha_t(ij)a
296			IF (first .AND. ini_anal) teta(ij, l) = a
297			END DO
298			END DO
299			END IF
300	guez	3
301	guez	29	! P
302			IF (guide_p) THEN
303			DO ij = 1, ip1jmp1
304			a = (1.-tau)psrea1(ij) + taupsrea2(ij)
305			ps(ij) = (1.-alpha_p(ij))ps(ij) + alpha_p(ij)a
306			IF (first .AND. ini_anal) ps(ij) = a
307			END DO
308	guez	37	forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
309	guez	29	CALL massdair(p, masse)
310			END IF
311	guez	3
312
313	guez	29	! q
314			IF (guide_q) THEN
315			DO l = 1, llm
316			DO ij = 1, ip1jmp1
317			a = (1.-tau)qrea1(ij, l) + tauqrea2(ij, l)
318			! hum relative en % -> hum specif
319			a = qsat(ij, l)a0.01
320			q(ij, l) = (1.-alpha_q(ij))q(ij, l) + alpha_q(ij)a
321			IF (first .AND. ini_anal) q(ij, l) = a
322			END DO
323			END DO
324			END IF
325	guez	3
326	guez	29	! vcov
327			IF (guide_v) THEN
328			DO l = 1, llm
329			DO ij = 1, ip1jm
330			a = (1.-tau)vcovrea1(ij, l) + tauvcovrea2(ij, l)
331			vcov(ij, l) = (1.-alpha_v(ij))vcov(ij, l) + alpha_v(ij)a
332			IF (first .AND. ini_anal) vcov(ij, l) = a
333			END DO
334			IF (first .AND. ini_anal) vcov(ij, l) = a
335			END DO
336			END IF
337	guez	3
338	guez	29	first = .FALSE.
339	guez	3
340	guez	29	END SUBROUTINE guide
341	guez	20
342			END MODULE guide_m