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