/[lmdze]/trunk/dyn3d/Guide/guide.f
ViewVC logotype

Diff of /trunk/dyn3d/Guide/guide.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.22  
changed lines
  Added in v.36

  ViewVC Help
Powered by ViewVC 1.1.21