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