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 comconst, ONLY : cpp, daysec, dtvr, kappa
    USE comvert, ONLY : ap, bp, preff, presnivs
    USE conf_gcm_m, ONLY : day_step, iperiod
    USE comgeom, ONLY : aire, rlatu, rlonv
    USE dimens_m, ONLY : jjm, llm
    USE exner_hyb_m, ONLY : exner_hyb
    USE inigrads_m, ONLY : inigrads
    use netcdf, only: nf90_nowrite, nf90_open, nf90_close
    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

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