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, &
         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

    ! 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

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