dyn3d/Guide/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 aire_min, aire_max

CONTAINS

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

    ! Author: F.Hourdin

    USE comconst, ONLY: cpp, daysec, dtvr, kappa
    USE comgeom, ONLY: aire, rlatu, rlonv
    USE conf_gcm_m, ONLY: day_step, iperiod
    use conf_guide_m, only: conf_guide, guide_u, guide_v, guide_t, guide_q, &
         guide_p, ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &
         tau_min_t, tau_max_t, tau_min_q, tau_max_q, tau_min_p, tau_max_p, &
         online
    USE dimens_m, ONLY: jjm, llm
    USE disvert_m, ONLY: ap, bp, preff, presnivs
    USE exner_hyb_m, ONLY: exner_hyb
    USE inigrads_m, ONLY: inigrads
    use massdair_m, only: massdair
    use netcdf, only: nf90_nowrite, nf90_open, nf90_close, nf90_inq_dimid, &
         nf90_inquire_dimension
    use nr_util, only: pi
    USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1, llmp1
    USE q_sat_m, ONLY: q_sat
    USE serre, ONLY: clat, clon
    use tau2alpha_m, only: tau2alpha, dxdys

    INTEGER, INTENT(IN):: itau

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

    ! Local:

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

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

    INTEGER ilon, ilat
    REAL factt, ztau(ip1jmp1)

    INTEGER ij, l
    INTEGER ncidpl, varidpl, status
    INTEGER rcod, rid
    REAL ditau, tau, a
    INTEGER, 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.
    CHARACTER(len=10) file
    INTEGER:: igrads = 2
    REAL:: dtgrads = 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)
    CALL exner_hyb(ps, p, pks, pk, pkf)
    tnat(:, :) = pk(:, :)*teta(:, :)/cpp
    unskap = 1./kappa
    pres(:, :) = preff*(pk(:, :)/cpp)**unskap
    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

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

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