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

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

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

revision 90 by guez, Wed Mar 12 21:16:36 2014 UTC revision 91 by guez, Wed Mar 26 17:18:58 2014 UTC
# Line 103  CONTAINS Line 103  CONTAINS
103      ! alpha=1 signifie pas d'injection      ! alpha=1 signifie pas d'injection
104      ! alpha=0 signifie injection totale      ! alpha=0 signifie injection totale
105    
106      IF (online== - 1) THEN      IF (online /= - 1) THEN
107         RETURN         IF (first) THEN
108      END IF            CALL conf_guide
109              file = 'guide'
110      IF (first) THEN            CALL inigrads(igrads, rlonv, 180. / pi, -180., 180., rlatu, -90., 90., &
111         CALL conf_guide                 180. / pi, presnivs, 1., dtgrads, file, 'dyn_zon ')
112         file = 'guide'            PRINT *, '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
113         CALL inigrads(igrads, rlonv, 180. / pi, -180., 180., rlatu, -90., 90., &  
114              180. / pi, presnivs, 1., dtgrads, file, 'dyn_zon ')            IF (online==1) THEN
115         PRINT *, '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'               ! Constantes de temps de rappel en jour
116         IF (online== - 1) RETURN               ! 0.1 c'est en gros 2h30.
117                 ! 1e10 est une constante infinie donc en gros pas de guidage
118         IF (online==1) THEN  
119            ! Constantes de temps de rappel en jour               ! coordonnees du centre du zoom
120            ! 0.1 c'est en gros 2h30.               CALL coordij(clon, clat, ilon, ilat)
121            ! 1e10 est une constante infinie donc en gros pas de guidage               ! aire de la maille au centre du zoom
122                 aire_min = aire(ilon+(ilat - 1) * iip1)
123            ! coordonnees du centre du zoom               ! aire maximale de la maille
124            CALL coordij(clon, clat, ilon, ilat)               aire_max = 0.
125            ! aire de la maille au centre du zoom               DO ij = 1, ip1jmp1
126            aire_min = aire(ilon+(ilat - 1) * iip1)                  aire_max = max(aire_max, aire(ij))
127            ! aire maximale de la maille               END DO
128            aire_max = 0.               ! factt = pas de temps en fraction de jour
129            DO ij = 1, ip1jmp1               factt = dtvr * iperiod / daysec
              aire_max = max(aire_max, aire(ij))  
           END DO  
           ! factt = pas de temps en fraction de jour  
           factt = dtvr * iperiod / daysec  
130    
131            CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)               CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
132            CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)               CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
133            CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)               CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)
134            CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)               CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)
135            CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)               CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)
136    
137            CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')               CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
138            CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')               CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')
139            CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')               CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')
140    
141                 ! Cas ou on force exactement par les variables analysees
142              ELSE
143                 alpha_t = 0.
144                 alpha_u = 0.
145                 alpha_v = 0.
146                 alpha_p = 0.
147                 ! physic=.false.
148              END IF
149    
150              itau_test = 1001
151              step_rea = 1
152              count_no_rea = 0
153              ncidpl = -99
154    
155              ! itau_test montre si l'importation a deja ete faite au rang itau
156              ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
157              if (guide_u) then
158                 if (ncidpl.eq. - 99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
159              endif
160    
161              if (guide_v) then
162                 if (ncidpl.eq. - 99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
163              endif
164    
165              if (guide_T) then
166                 if (ncidpl.eq. - 99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
167              endif
168    
169              if (guide_Q) then
170                 if (ncidpl.eq. - 99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
171              endif
172    
173              IF (ncep) THEN
174                 status = nf90_inq_dimid(ncidpl, 'LEVEL', rid)
175              ELSE
176                 status = nf90_inq_dimid(ncidpl, 'PRESSURE', rid)
177              END IF
178              status = nf90_inquire_dimension(ncidpl, rid, len=nlev)
179              PRINT *, 'nlev', nlev
180              rcod = nf90_close(ncidpl)
181              ! Lecture du premier etat des reanalyses.
182              CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
183                   masserea2, nlev)
184              qrea2 = max(qrea2, 0.1)
185    
186            ! Cas ou on force exactement par les variables analysees            ! Debut de l'integration temporelle:
187         ELSE         END IF ! first
           alpha_t = 0.  
           alpha_u = 0.  
           alpha_v = 0.  
           alpha_p = 0.  
           ! physic=.false.  
        END IF  
188    
189         itau_test = 1001         ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
        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  
190    
191         IF (ncep) THEN         ditau = real(itau)
192            status = nf90_inq_dimid(ncidpl, 'LEVEL', rid)         dday_step = real(day_step)
193           WRITE (*, *) 'ditau, dday_step'
194           WRITE (*, *) ditau, dday_step
195           toto = 4 * ditau / dday_step
196           reste = toto - aint(toto)
197    
198           IF (reste==0.) THEN
199              IF (itau_test==itau) THEN
200                 WRITE (*, *) 'deuxieme passage de advreel a itau=', itau
201                 STOP
202              ELSE
203                 vcovrea1 = vcovrea2
204                 ucovrea1 = ucovrea2
205                 tetarea1 = tetarea2
206                 qrea1 = qrea2
207    
208                 PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
209                      count_no_rea, ' non lectures'
210                 step_rea = step_rea + 1
211                 itau_test = itau
212                 CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &
213                      qrea2, masserea2, nlev)
214                 qrea2 = max(qrea2, 0.1)
215                 factt = dtvr * iperiod / daysec
216                 ztau = factt / max(alpha_t, 1E-10)
217                 CALL wrgrads(igrads, 1, aire, 'aire ', 'aire ')
218                 CALL wrgrads(igrads, 1, dxdys, 'dxdy ', 'dxdy ')
219                 CALL wrgrads(igrads, 1, alpha_u, 'au ', 'au ')
220                 CALL wrgrads(igrads, 1, alpha_t, 'at ', 'at ')
221                 CALL wrgrads(igrads, 1, ztau, 'taut ', 'taut ')
222                 CALL wrgrads(igrads, llm, ucov, 'u ', 'u ')
223                 CALL wrgrads(igrads, llm, ucovrea2, 'ua ', 'ua ')
224                 CALL wrgrads(igrads, llm, teta, 'T ', 'T ')
225                 CALL wrgrads(igrads, llm, tetarea2, 'Ta ', 'Ta ')
226                 CALL wrgrads(igrads, llm, qrea2, 'Qa ', 'Qa ')
227                 CALL wrgrads(igrads, llm, q, 'Q ', 'Q ')
228    
229                 CALL wrgrads(igrads, llm, qsat, 'QSAT ', 'QSAT ')
230    
231              END IF
232         ELSE         ELSE
233            status = nf90_inq_dimid(ncidpl, 'PRESSURE', rid)            count_no_rea = count_no_rea + 1
234         END IF         END IF
        status = nf90_inquire_dimension(ncidpl, rid, len=nlev)  
        PRINT *, 'nlev', nlev  
        rcod = nf90_close(ncidpl)  
        ! Lecture du premier etat des reanalyses.  
        CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &  
             masserea2, nlev)  
        qrea2 = max(qrea2, 0.1)  
   
        ! Debut de l'integration temporelle:  
     END IF ! first  
   
     ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:  
   
     ditau = real(itau)  
     dday_step = real(day_step)  
     WRITE (*, *) 'ditau, dday_step'  
     WRITE (*, *) ditau, dday_step  
     toto = 4 * ditau / dday_step  
     reste = toto - aint(toto)  
   
     IF (reste==0.) THEN  
        IF (itau_test==itau) THEN  
           WRITE (*, *) 'deuxieme passage de advreel a itau=', itau  
           STOP  
        ELSE  
           vcovrea1 = vcovrea2  
           ucovrea1 = ucovrea2  
           tetarea1 = tetarea2  
           qrea1 = qrea2  
   
           PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &  
                count_no_rea, ' non lectures'  
           step_rea = step_rea + 1  
           itau_test = itau  
           CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &  
                qrea2, masserea2, nlev)  
           qrea2 = max(qrea2, 0.1)  
           factt = dtvr * iperiod / daysec  
           ztau = factt / max(alpha_t, 1E-10)  
           CALL wrgrads(igrads, 1, aire, 'aire ', 'aire ')  
           CALL wrgrads(igrads, 1, dxdys, 'dxdy ', 'dxdy ')  
           CALL wrgrads(igrads, 1, alpha_u, 'au ', 'au ')  
           CALL wrgrads(igrads, 1, alpha_t, 'at ', 'at ')  
           CALL wrgrads(igrads, 1, ztau, 'taut ', 'taut ')  
           CALL wrgrads(igrads, llm, ucov, 'u ', 'u ')  
           CALL wrgrads(igrads, llm, ucovrea2, 'ua ', 'ua ')  
           CALL wrgrads(igrads, llm, teta, 'T ', 'T ')  
           CALL wrgrads(igrads, llm, tetarea2, 'Ta ', 'Ta ')  
           CALL wrgrads(igrads, llm, qrea2, 'Qa ', 'Qa ')  
           CALL wrgrads(igrads, llm, q, 'Q ', 'Q ')  
235    
236            CALL wrgrads(igrads, llm, qsat, 'QSAT ', 'QSAT ')         ! Guidage
237           ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
238    
239           IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'
240    
241           ditau = real(itau)
242           dday_step = real(day_step)
243    
244           tau = 4 * ditau / dday_step
245           tau = tau - aint(tau)
246    
247           ! ucov
248           IF (guide_u) THEN
249              DO l = 1, llm
250                 DO ij = 1, ip1jmp1
251                    a = (1. - tau) * ucovrea1(ij, l) + tau * ucovrea2(ij, l)
252                    ucov(ij, l) = (1. - alpha_u(ij)) * ucov(ij, l) + alpha_u(ij) * a
253                    IF (first .AND. ini_anal) ucov(ij, l) = a
254                 END DO
255              END DO
256           END IF
257    
258           IF (guide_t) THEN
259              DO l = 1, llm
260                 do j = 1, jjm + 1
261                    DO i = 1, iim + 1
262                       a = (1. - tau) * tetarea1(i, j, l) + tau * tetarea2(i, j, l)
263                       teta(i, j, l) = (1. - alpha_t(i, j)) * teta(i, j, l) &
264                            + alpha_t(i, j) * a
265                       IF (first .AND. ini_anal) teta(i, j, l) = a
266                    END DO
267                 end do
268              END DO
269         END IF         END IF
270      ELSE  
271         count_no_rea = count_no_rea + 1         IF (guide_q) THEN
272      END IF            DO l = 1, llm
273                 do j = 1, jjm + 1
274      ! Guidage                  DO i = 1, iim + 1
275      ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses                     a = (1. - tau) * qrea1(i, j, l) + tau * qrea2(i, j, l)
276                       ! hum relative en % -> hum specif
277      IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'                     a = qsat(i, j, l) * a * 0.01
278                       q(i, j, l) = (1. - alpha_q(i, j)) * q(i, j, l) &
279      ditau = real(itau)                          + alpha_q(i, j) * a
280      dday_step = real(day_step)                     IF (first .AND. ini_anal) q(i, j, l) = a
281                    END DO
282      tau = 4 * ditau / dday_step               end do
     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  
283            END DO            END DO
284         END DO         END IF
     END IF  
285    
286      IF (guide_t) THEN         ! vcov
287         DO l = 1, llm         IF (guide_v) THEN
288            do j = 1, jjm + 1            DO l = 1, llm
289               DO i = 1, iim + 1               DO ij = 1, ip1jm
290                  a = (1. - tau) * tetarea1(i, j, l) + tau * tetarea2(i, j, l)                  a = (1. - tau) * vcovrea1(ij, l) + tau * vcovrea2(ij, l)
291                  teta(i, j, l) = (1. - alpha_t(i, j)) * teta(i, j, l) &                  vcov(ij, l) = (1. - alpha_v(ij)) * vcov(ij, l) + alpha_v(ij) * a
292                       + alpha_t(i, j) * a                  IF (first .AND. ini_anal) vcov(ij, l) = a
                 IF (first .AND. ini_anal) teta(i, j, l) = a  
293               END DO               END DO
           end do  
        END DO  
     END IF  
   
     IF (guide_q) THEN  
        DO l = 1, llm  
           do j = 1, jjm + 1  
              DO i = 1, iim + 1  
                 a = (1. - tau) * qrea1(i, j, l) + tau * qrea2(i, j, l)  
                 ! hum relative en % -> hum specif  
                 a = qsat(i, j, l) * a * 0.01  
                 q(i, j, l) = (1. - alpha_q(i, j)) * q(i, j, l) &  
                      + alpha_q(i, j) * a  
                 IF (first .AND. ini_anal) q(i, j, l) = a  
              END DO  
           end do  
        END DO  
     END IF  
   
     ! vcov  
     IF (guide_v) THEN  
        DO l = 1, llm  
           DO ij = 1, ip1jm  
              a = (1. - tau) * vcovrea1(ij, l) + tau * vcovrea2(ij, l)  
              vcov(ij, l) = (1. - alpha_v(ij)) * vcov(ij, l) + alpha_v(ij) * a  
294               IF (first .AND. ini_anal) vcov(ij, l) = a               IF (first .AND. ini_anal) vcov(ij, l) = a
295            END DO            END DO
296            IF (first .AND. ini_anal) vcov(ij, l) = a         END IF
        END DO  
     END IF  
297    
298      first = .FALSE.         first = .FALSE.
299        end IF
300    
301    END SUBROUTINE guide    END SUBROUTINE guide
302    
Legend:
Removed from v.90  
changed lines
  Added in v.91
  ViewVC Help
Powered by ViewVC 1.1.21