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 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, &
         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: iim, 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 read_reanalyse_m, only: read_reanalyse
    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(iim + 1, jjm + 1, llm) ! température potentielle 
    REAL q(iim + 1, jjm + 1, llm)
    REAL, intent(out):: masse(ip1jmp1, llm) ! masse d'air
    REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol

    ! Local:

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

    REAL, save:: alpha_q(iim + 1, jjm + 1)
    REAL, save:: alpha_t(iim + 1, jjm + 1), 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(iim + 1, jjm + 1)

    INTEGER ij, i, j, l
    INTEGER ncidpl, status
    INTEGER rcod, rid
    REAL ditau, tau, a
    INTEGER, SAVE:: nlev

    ! TEST SUR QSAT
    REAL p(iim + 1, jjm + 1, llmp1)
    real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
    REAL pres(iim + 1, jjm + 1, llm)

    REAL qsat(iim + 1, jjm + 1, llm)
    REAL unskap
    REAL tnat(iim + 1, jjm + 1, 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)
    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
       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) 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, 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, nlev)
             qrea2 = max(qrea2, 0.1)
             factt = dtvr * iperiod / daysec
             ztau = factt / max(alpha_t, 1E-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 j = 1, jjm + 1
                DO i = 1, iim + 1
                   a = (1. - tau) * tetarea1(i, j, l) + tau * tetarea2(i, j, l)
                   teta(i, j, l) = (1. - alpha_t(i, j)) * teta(i, j, l) &
                        + alpha_t(i, j) * a
                   IF (first .AND. ini_anal) teta(i, j, l) = a
                END DO
             end do
          END DO
       END IF

       IF (guide_q) THEN
          DO l = 1, llm
             do j = 1, jjm + 1
                DO i = 1, iim + 1
                   a = (1. - tau) * qrea1(i, j, l) + tau * qrea2(i, j, l)
                   ! hum relative en % -> hum specif
                   a = qsat(i, j, l) * a * 0.01
                   q(i, j, l) = (1. - alpha_q(i, j)) * q(i, j, l) &
                        + alpha_q(i, j) * a
                   IF (first .AND. ini_anal) q(i, j, l) = a
                END DO
             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 IF

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