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

Diff of /trunk/dyn3d/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 29 by guez, Tue Mar 30 10:44:42 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
# Line 19  MODULE guide_m Line 19  MODULE guide_m
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        IMPLICIT NONE
40    
41        INCLUDE 'netcdf.inc'
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        !   common passe pour des sorties
51        REAL :: dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
52        COMMON /comdxdy/dxdys, dxdyu, dxdyv
53    
54        !   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    
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           PRINT *, 'initialisation du guide '
138           CALL conf_guide
139           PRINT *, 'apres conf_guide'
140    
141           file = 'guide'
142           CALL inigrads(igrads, rlonv, 180./pi, -180., 180., rlatu, -90., 90., &
143                180./pi, presnivs, 1., dtgrads, file, 'dyn_zon ')
144    
145           PRINT *, &
146                '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
147    
148           IF (online==-1) RETURN
149           IF (online==1) THEN
150    
151              !  Constantes de temps de rappel en jour
152              !  0.1 c'est en gros 2h30.
153              !  1e10  est une constante infinie donc en gros pas de guidage
154    
155              !   coordonnees du centre du zoom
156              CALL coordij(clon, clat, ilon, ilat)
157              !   aire de la maille au centre du zoom
158              aire_min = aire(ilon+(ilat-1)*iip1)
159              !   aire maximale de la maille
160              aire_max = 0.
161              DO ij = 1, ip1jmp1
162                 aire_max = max(aire_max, aire(ij))
163              END DO
164              !  factt = pas de temps en fraction de jour
165              factt = dtvr*iperiod/daysec
166    
167              CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
168              CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
169              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              CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)
172    
173              CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
174              CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U   ')
175              CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T   ')
176    
177              !   Cas ou on force exactement par les variables analysees
178           ELSE
179              alpha_t = 0.
180              alpha_u = 0.
181              alpha_v = 0.
182              alpha_p = 0.
183              !           physic=.false.
184           END IF
185    
186           itau_test = 1001
187           step_rea = 1
188           count_no_rea = 0
189           ncidpl = -99
190    
191           !    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           if (guide_u) then
194              if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
195           endif
196    
197           if (guide_v) then
198              if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
199           endif
200    
201           if (guide_T) then
202              if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
203           endif
204    
205           if (guide_Q) then
206              if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
207           endif
208    
209           IF (ncep) THEN
210              status = nf_inq_dimid(ncidpl, 'LEVEL', rid)
211           ELSE
212              status = nf_inq_dimid(ncidpl, 'PRESSURE', rid)
213           END IF
214           status = nf_inq_dimlen(ncidpl, rid, nlev)
215           PRINT *, 'nlev', nlev
216           rcod = nf90_close(ncidpl)
217           !   Lecture du premier etat des reanalyses.
218           CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
219                masserea2, psrea2, 1, nlev)
220           qrea2(:, :) = max(qrea2(:, :), 0.1)
221    
222    
223           !   Debut de l'integration temporelle:
224        END IF ! first
225    
226        ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
227    
228        ditau = real(itau)
229        dday_step = real(day_step)
230        WRITE (*, *) 'ditau, dday_step'
231        WRITE (*, *) ditau, dday_step
232        toto = 4*ditau/dday_step
233        reste = toto - aint(toto)
234    
235        IF (reste==0.) THEN
236           IF (itau_test==itau) THEN
237              WRITE (*, *) 'deuxieme passage de advreel a itau=', itau
238              STOP
239           ELSE
240              vcovrea1(:, :) = vcovrea2(:, :)
241              ucovrea1(:, :) = ucovrea2(:, :)
242              tetarea1(:, :) = tetarea2(:, :)
243              qrea1(:, :) = qrea2(:, :)
244    
245              PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
246                   count_no_rea, ' non lectures'
247              step_rea = step_rea + 1
248              itau_test = itau
249              CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &
250                   qrea2, masserea2, psrea2, 1, nlev)
251              qrea2(:, :) = max(qrea2(:, :), 0.1)
252              factt = dtvr*iperiod/daysec
253              ztau(:) = factt/max(alpha_t(:), 1.E-10)
254              CALL wrgrads(igrads, 1, aire, 'aire      ', 'aire      ')
255              CALL wrgrads(igrads, 1, dxdys, 'dxdy      ', 'dxdy      ')
256              CALL wrgrads(igrads, 1, alpha_u, 'au        ', 'au        ')
257              CALL wrgrads(igrads, 1, alpha_t, 'at        ', 'at        ')
258              CALL wrgrads(igrads, 1, ztau, 'taut      ', 'taut      ')
259              CALL wrgrads(igrads, llm, ucov, 'u         ', 'u         ')
260              CALL wrgrads(igrads, llm, ucovrea2, 'ua        ', 'ua        ')
261              CALL wrgrads(igrads, llm, teta, 'T         ', 'T         ')
262              CALL wrgrads(igrads, llm, tetarea2, 'Ta        ', 'Ta        ')
263              CALL wrgrads(igrads, llm, qrea2, 'Qa        ', 'Qa        ')
264              CALL wrgrads(igrads, llm, q, 'Q         ', 'Q         ')
265    
266              CALL wrgrads(igrads, llm, qsat, 'QSAT      ', 'QSAT      ')
267    
268           END IF
269        ELSE
270           count_no_rea = count_no_rea + 1
271        END IF
272    
273        !   Guidage
274        !    x_gcm = a * x_gcm + (1-a) * x_reanalyses
275    
276        IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'
277    
278        ditau = real(itau)
279        dday_step = real(day_step)
280    
281    
282        tau = 4*ditau/dday_step
283        tau = tau - aint(tau)
284    
285        !  ucov
286        IF (guide_u) THEN
287           DO l = 1, llm
288              DO ij = 1, ip1jmp1
289                 a = (1.-tau)*ucovrea1(ij, l) + tau*ucovrea2(ij, l)
290                 ucov(ij, l) = (1.-alpha_u(ij))*ucov(ij, l) + alpha_u(ij)*a
291                 IF (first .AND. ini_anal) ucov(ij, l) = a
292              END DO
293           END DO
294        END IF
295    
296        IF (guide_t) THEN
297           DO l = 1, llm
298              DO ij = 1, ip1jmp1
299                 a = (1.-tau)*tetarea1(ij, l) + tau*tetarea2(ij, l)
300                 teta(ij, l) = (1.-alpha_t(ij))*teta(ij, l) + alpha_t(ij)*a
301                 IF (first .AND. ini_anal) teta(ij, l) = a
302              END DO
303           END DO
304        END IF
305    
306        !  P
307        IF (guide_p) THEN
308           DO ij = 1, ip1jmp1
309              a = (1.-tau)*psrea1(ij) + tau*psrea2(ij)
310              ps(ij) = (1.-alpha_p(ij))*ps(ij) + alpha_p(ij)*a
311              IF (first .AND. ini_anal) ps(ij) = a
312           END DO
313           CALL pression(ip1jmp1, ap, bp, ps, p)
314           CALL massdair(p, masse)
315        END IF
316    
317    
318        !  q
319        IF (guide_q) THEN
320           DO l = 1, llm
321              DO ij = 1, ip1jmp1
322                 a = (1.-tau)*qrea1(ij, l) + tau*qrea2(ij, l)
323                 !   hum relative en % -> hum specif
324                 a = qsat(ij, l)*a*0.01
325                 q(ij, l) = (1.-alpha_q(ij))*q(ij, l) + alpha_q(ij)*a
326                 IF (first .AND. ini_anal) q(ij, l) = a
327              END DO
328           END DO
329        END IF
330    
331        ! vcov
332        IF (guide_v) THEN
333           DO l = 1, llm
334              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              END DO
339              IF (first .AND. ini_anal) vcov(ij, l) = a
340           END DO
341        END IF
342    
343    !      ......   Version  du 10/01/98    ..........      first = .FALSE.
344    
345    !             avec  coordonnees  verticales hybrides    END SUBROUTINE guide
   !   avec nouveaux operat. dissipation * ( gradiv2, divgrad2, nxgraro2 )  
   
   !=======================================================================  
   
   !   Auteur:  F.Hourdin  
   !   -------  
   
   !   Objet:  
   !   ------  
   
   !   GCM LMD nouvelle grille  
   
   !=======================================================================  
   
   ! Dans inigeom , nouveaux calculs pour les elongations  cu , cv  
   ! et possibilite d'appeler une fonction f(y)  a derivee tangente  
   ! hyperbolique a la  place de la fonction a derivee sinusoidale.          
   
   !   ...  Possibilite de choisir le shema de Van-leer pour l'advection de  
   !         q  , en faisant iadv = 10  dans   traceur  (29/04/97) .  
   
   !-----------------------------------------------------------------------  
   !   Declarations:  
   !   -------------  
   
   
   !   variables dynamiques  
   REAL :: vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants  
   REAL :: 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  
   
   !IM 180305   real aire_min, aire_max  
   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)  
   !cccccccccccccccc  
   
   
   LOGICAL :: first  
   SAVE first  
   DATA 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  
   !-----------------------------------------------------------------------  
   CALL pression(ip1jmp1, ap, bp, ps, p)  
   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  
   
         !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
         !  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  
         !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
         !   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  
   
         !     subroutine tau2alpha(type, im, jm, factt, taumin, taumax, alpha)  
         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   ')  
   
         !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
         !   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)  
   !     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           ')  
   
   first = .FALSE.  
   
   RETURN  
 END SUBROUTINE guide  
346    
347    !=======================================================================    !=======================================================================
348    SUBROUTINE tau2alpha(type, pim, pjm, factt, taumin, taumax, alpha)    SUBROUTINE tau2alpha(type, pim, pjm, factt, taumin, taumax, alpha)

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

  ViewVC Help
Powered by ViewVC 1.1.21