/[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 91 by guez, Wed Mar 26 17:18:58 2014 UTC revision 107 by guez, Thu Sep 11 15:09:15 2014 UTC
# Line 9  MODULE guide_m Line 9  MODULE guide_m
9    
10  CONTAINS  CONTAINS
11    
12    SUBROUTINE guide(itau, ucov, vcov, teta, q, masse, ps)    SUBROUTINE guide(itau, ucov, vcov, teta, q, ps)
13    
14      ! Author: F.Hourdin      ! Author: F.Hourdin
15    
# Line 18  CONTAINS Line 18  CONTAINS
18      USE conf_gcm_m, ONLY: day_step, iperiod      USE conf_gcm_m, ONLY: day_step, iperiod
19      use conf_guide_m, only: conf_guide, guide_u, guide_v, guide_t, guide_q, &      use conf_guide_m, only: conf_guide, guide_u, guide_v, guide_t, guide_q, &
20           ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &           ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &
21           tau_min_t, tau_max_t, tau_min_q, tau_max_q, tau_min_p, tau_max_p, &           tau_min_t, tau_max_t, tau_min_q, tau_max_q, online
          online  
22      USE dimens_m, ONLY: iim, jjm, llm      USE dimens_m, ONLY: iim, jjm, llm
23      USE disvert_m, ONLY: ap, bp, preff, presnivs      USE disvert_m, ONLY: ap, bp, preff, presnivs
24        use dump2d_m, only: dump2d
25      USE exner_hyb_m, ONLY: exner_hyb      USE exner_hyb_m, ONLY: exner_hyb
26      USE inigrads_m, ONLY: inigrads      USE inigrads_m, ONLY: inigrads
27      use massdair_m, only: massdair      use netcdf, only: nf90_nowrite
28      use netcdf, only: nf90_nowrite, nf90_open, nf90_close, nf90_inq_dimid, &      use netcdf95, only: nf95_close, nf95_inq_dimid, nf95_inquire_dimension, &
29           nf90_inquire_dimension           nf95_open
30      use nr_util, only: pi      use nr_util, only: pi
31      USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1, llmp1      USE paramet_m, ONLY: iip1, ip1jmp1, jjp1, llmp1
32      USE q_sat_m, ONLY: q_sat      USE q_sat_m, ONLY: q_sat
33      use read_reanalyse_m, only: read_reanalyse      use read_reanalyse_m, only: read_reanalyse
34      USE serre, ONLY: clat, clon      USE serre, ONLY: clat, clon
# Line 37  CONTAINS Line 37  CONTAINS
37      INTEGER, INTENT(IN):: itau      INTEGER, INTENT(IN):: itau
38    
39      ! variables dynamiques      ! variables dynamiques
40      REAL ucov(ip1jmp1, llm), vcov(ip1jm, llm) ! vents covariants  
41        REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
42        REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
43    
44      REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! température potentielle      REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! température potentielle
45      REAL q(iim + 1, jjm + 1, llm)      REAL, intent(inout):: q(iim + 1, jjm + 1, llm)
     REAL, intent(out):: masse(ip1jmp1, llm) ! masse d'air  
46      REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol      REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol
47    
48      ! Local:      ! Local:
49    
50      ! variables dynamiques pour les reanalyses.      ! variables dynamiques pour les reanalyses.
51      REAL, save:: ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas  
52        REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
53        ! vents covariants reanalyses
54    
55      REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales      REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales
56      REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales      REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales
57      REAL, save:: ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas  
58        REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
59        ! vents covariants reanalyses
60    
61      REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales      REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales
62      REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales      REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales
63      REAL, save:: masserea2(ip1jmp1, llm) ! masse      REAL, save:: masserea2(ip1jmp1, llm) ! masse
64    
65        ! alpha determine la part des injections de donnees a chaque etape
66        ! alpha=1 signifie pas d'injection
67        ! alpha=0 signifie injection totale
68      REAL, save:: alpha_q(iim + 1, jjm + 1)      REAL, save:: alpha_q(iim + 1, jjm + 1)
69      REAL, save:: alpha_t(iim + 1, jjm + 1), alpha_p(ip1jmp1)      REAL, save:: alpha_t(iim + 1, jjm + 1)
70      REAL, save:: alpha_u(ip1jmp1), alpha_v(ip1jm)      REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)
71      REAL dday_step, toto, reste  
     real, save:: itau_test  
72      INTEGER, save:: step_rea, count_no_rea      INTEGER, save:: step_rea, count_no_rea
73    
74      INTEGER ilon, ilat      INTEGER ilon, ilat
75      REAL factt, ztau(iim + 1, jjm + 1)      REAL factt ! pas de temps entre deux appels au guidage, en fraction de jour
76        real ztau(iim + 1, jjm + 1)
77    
78      INTEGER ij, i, j, l      INTEGER ij, l
79      INTEGER ncidpl, status      INTEGER ncidpl
80      INTEGER rcod, rid      INTEGER rid
81      REAL ditau, tau, a      REAL tau
82      INTEGER, SAVE:: nlev      INTEGER, SAVE:: nlev
83    
84      ! TEST SUR QSAT      ! TEST SUR QSAT
85      REAL p(iim + 1, jjm + 1, llmp1)      REAL p(iim + 1, jjm + 1, llmp1)
86      real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)      real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
     REAL pres(iim + 1, jjm + 1, llm)  
87    
88      REAL qsat(iim + 1, jjm + 1, llm)      REAL qsat(iim + 1, jjm + 1, llm)
     REAL unskap  
     REAL tnat(iim + 1, jjm + 1, llm)  
89    
90      LOGICAL:: first = .TRUE.      INTEGER, parameter:: igrads = 2
     CHARACTER(len=10) file  
     INTEGER:: igrads = 2  
91      REAL:: dtgrads = 100.      REAL:: dtgrads = 100.
92    
93      !-----------------------------------------------------------------------      !-----------------------------------------------------------------------
94    
95      PRINT *, 'Call sequence information: guide'      PRINT *, 'Call sequence information: guide'
96    
97      ! calcul de l'humidite saturante      first_call: IF (itau == 0) THEN
98           CALL conf_guide
99      forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps         CALL inigrads(igrads, rlonv, 180. / pi, -180., 180., rlatu, -90., &
100      CALL massdair(p, masse)              90., 180. / pi, presnivs, 1., dtgrads, 'guide', 'dyn_zon ')
101      CALL exner_hyb(ps, p, pks, pk)  
102      tnat = pk * teta / cpp         IF (online) THEN
103      unskap = 1. / kappa            ! Constantes de temps de rappel en jour
104      pres = preff * (pk / cpp)**unskap  
105      qsat = q_sat(tnat, pres)            ! coordonnees du centre du zoom
106              CALL coordij(clon, clat, ilon, ilat)
107      ! initialisations pour la lecture des reanalyses.            ! aire de la maille au centre du zoom
108      ! alpha determine la part des injections de donnees a chaque etape            aire_min = aire(ilon+(ilat - 1) * iip1)
109      ! alpha=1 signifie pas d'injection            ! aire maximale de la maille
110      ! alpha=0 signifie injection totale            aire_max = 0.
111              DO ij = 1, ip1jmp1
112      IF (online /= - 1) THEN               aire_max = max(aire_max, aire(ij))
113         IF (first) THEN            END DO
           CALL conf_guide  
           file = 'guide'  
           CALL inigrads(igrads, rlonv, 180. / pi, -180., 180., rlatu, -90., 90., &  
                180. / pi, presnivs, 1., dtgrads, file, 'dyn_zon ')  
           PRINT *, '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'  
   
           IF (online==1) THEN  
              ! Constantes de temps de rappel en jour  
              ! 0.1 c'est en gros 2h30.  
              ! 1e10 est une constante infinie donc en gros pas de guidage  
   
              ! coordonnees du centre du zoom  
              CALL coordij(clon, clat, ilon, ilat)  
              ! aire de la maille au centre du zoom  
              aire_min = aire(ilon+(ilat - 1) * iip1)  
              ! aire maximale de la maille  
              aire_max = 0.  
              DO ij = 1, ip1jmp1  
                 aire_max = max(aire_max, aire(ij))  
              END DO  
              ! factt = pas de temps en fraction de jour  
              factt = dtvr * iperiod / daysec  
   
              CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)  
              CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)  
              CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)  
              CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)  
              CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)  
   
              CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')  
              CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')  
              CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')  
   
              ! Cas ou on force exactement par les variables analysees  
           ELSE  
              alpha_t = 0.  
              alpha_u = 0.  
              alpha_v = 0.  
              alpha_p = 0.  
              ! physic=.false.  
           END IF  
   
           itau_test = 1001  
           step_rea = 1  
           count_no_rea = 0  
           ncidpl = -99  
   
           ! itau_test montre si l'importation a deja ete faite au rang itau  
           ! lecture d'un fichier netcdf pour determiner le nombre de niveaux  
           if (guide_u) then  
              if (ncidpl.eq. - 99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)  
           endif  
   
           if (guide_v) then  
              if (ncidpl.eq. - 99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)  
           endif  
   
           if (guide_T) then  
              if (ncidpl.eq. - 99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)  
           endif  
   
           if (guide_Q) then  
              if (ncidpl.eq. - 99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)  
           endif  
   
           IF (ncep) THEN  
              status = nf90_inq_dimid(ncidpl, 'LEVEL', rid)  
           ELSE  
              status = nf90_inq_dimid(ncidpl, 'PRESSURE', rid)  
           END IF  
           status = nf90_inquire_dimension(ncidpl, rid, len=nlev)  
           PRINT *, 'nlev', nlev  
           rcod = nf90_close(ncidpl)  
           ! Lecture du premier etat des reanalyses.  
           CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &  
                masserea2, 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 ')  
114    
115               CALL wrgrads(igrads, llm, qsat, 'QSAT ', 'QSAT ')            factt = dtvr * iperiod / daysec
116    
117            END IF            CALL tau2alpha(3, factt, tau_min_v, tau_max_v, alpha_v)
118              CALL tau2alpha(2, factt, tau_min_u, tau_max_u, alpha_u)
119              CALL tau2alpha(1, factt, tau_min_t, tau_max_t, alpha_t)
120              CALL tau2alpha(1, factt, tau_min_q, tau_max_q, alpha_q)
121    
122              CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
123              CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')
124              CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')
125         ELSE         ELSE
126            count_no_rea = count_no_rea + 1            ! Cas où on force exactement par les variables analysées
127              alpha_t = 0.
128              alpha_u = 0.
129              alpha_v = 0.
130              alpha_q = 0.
131         END IF         END IF
132    
133         ! Guidage         step_rea = 1
134         ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses         count_no_rea = 0
135           ncidpl = -99
136         IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'  
137           ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
138         ditau = real(itau)         if (guide_u) call nf95_open('u.nc',Nf90_NOWRITe,ncidpl)
139         dday_step = real(day_step)         if (guide_v) call nf95_open('v.nc',nf90_nowrite,ncidpl)
140           if (guide_T) call nf95_open('T.nc',nf90_nowrite,ncidpl)
141         tau = 4 * ditau / dday_step         if (guide_Q) call nf95_open('hur.nc',nf90_nowrite, ncidpl)
        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  
142    
143         IF (guide_t) THEN         IF (ncep) THEN
144            DO l = 1, llm            call nf95_inq_dimid(ncidpl, 'LEVEL', rid)
145               do j = 1, jjm + 1         ELSE
146                  DO i = 1, iim + 1            call nf95_inq_dimid(ncidpl, 'PRESSURE', rid)
                    a = (1. - tau) * tetarea1(i, j, l) + tau * tetarea2(i, j, l)  
                    teta(i, j, l) = (1. - alpha_t(i, j)) * teta(i, j, l) &  
                         + alpha_t(i, j) * a  
                    IF (first .AND. ini_anal) teta(i, j, l) = a  
                 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  
                 IF (first .AND. ini_anal) vcov(ij, l) = a  
              END DO  
              IF (first .AND. ini_anal) vcov(ij, l) = a  
           END DO  
147         END IF         END IF
148           call nf95_inquire_dimension(ncidpl, rid, nclen=nlev)
149         first = .FALSE.         PRINT *, 'nlev', nlev
150      end IF         call nf95_close(ncidpl)
151           ! Lecture du premier etat des reanalyses.
152           CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
153                masserea2, nlev)
154           qrea2 = max(qrea2, 0.1)
155        END IF first_call
156    
157        ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
158    
159        ! Nudging fields are given 4 times per day:
160        IF (mod(itau, day_step / 4) == 0) THEN
161           vcovrea1 = vcovrea2
162           ucovrea1 = ucovrea2
163           tetarea1 = tetarea2
164           qrea1 = qrea2
165    
166           PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
167                count_no_rea, ' non lectures'
168           step_rea = step_rea + 1
169           CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
170                masserea2, nlev)
171           qrea2 = max(qrea2, 0.1)
172           factt = dtvr * iperiod / daysec
173           ztau = factt / max(alpha_t, 1E-10)
174           CALL wrgrads(igrads, 1, aire, 'aire ', 'aire ')
175           CALL wrgrads(igrads, 1, dxdys, 'dxdy ', 'dxdy ')
176           CALL wrgrads(igrads, 1, alpha_u, 'au ', 'au ')
177           CALL wrgrads(igrads, 1, alpha_t, 'at ', 'at ')
178           CALL wrgrads(igrads, 1, ztau, 'taut ', 'taut ')
179           CALL wrgrads(igrads, llm, ucov, 'u ', 'u ')
180           CALL wrgrads(igrads, llm, ucovrea2, 'ua ', 'ua ')
181           CALL wrgrads(igrads, llm, teta, 'T ', 'T ')
182           CALL wrgrads(igrads, llm, tetarea2, 'Ta ', 'Ta ')
183           CALL wrgrads(igrads, llm, qrea2, 'Qa ', 'Qa ')
184           CALL wrgrads(igrads, llm, q, 'Q ', 'Q ')
185        ELSE
186           count_no_rea = count_no_rea + 1
187        END IF
188    
189        ! Guidage
190    
191        tau = mod(real(itau) / real(day_step / 4), 1.)
192    
193        ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
194    
195        IF (guide_u) THEN
196           IF (itau == 0 .AND. ini_anal) then
197              ucov = ucovrea1
198           else
199              forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) * ucov(:, :, l) &
200                   + alpha_u * ((1. - tau) * ucovrea1(:, :, l) &
201                   + tau * ucovrea2(:, :, l))
202           end IF
203        END IF
204    
205        IF (guide_t) THEN
206           IF (itau == 0 .AND. ini_anal) then
207              teta = tetarea1
208           else
209              forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) * teta(:, :, l) &
210                   + alpha_t * ((1. - tau) * tetarea1(:, :, l) &
211                   + tau * tetarea2(:, :, l))
212           end IF
213        END IF
214    
215        IF (guide_q) THEN
216           ! Calcul de l'humidité saturante :
217           forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
218           CALL exner_hyb(ps, p, pks, pk)
219           qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
220    
221           ! humidité relative en % -> humidité spécifique
222           IF (itau == 0 .AND. ini_anal) then
223              q = qsat * qrea1 * 0.01
224           else
225              forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
226                   + alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
227                   + tau * qrea2(:, :, l)) * 0.01)
228           end IF
229        END IF
230    
231        IF (guide_v) THEN
232           IF (itau == 0 .AND. ini_anal) then
233              vcov = vcovrea1
234           else
235              forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) * vcov(:, :, l) &
236                   + alpha_v * ((1. - tau) * vcovrea1(:, :, l) &
237                   + tau * vcovrea2(:, :, l))
238           end IF
239        END IF
240    
241    END SUBROUTINE guide    END SUBROUTINE guide
242    
Legend:
Removed from v.91  
changed lines
  Added in v.107
  ViewVC Help
Powered by ViewVC 1.1.21