/[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 109 by guez, Wed Sep 17 10:08:00 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
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    
75      INTEGER ij, i, j, l      INTEGER ij, l
76      INTEGER ncidpl, status      INTEGER ncid, dimid
77      INTEGER rcod, rid      REAL tau
     REAL ditau, tau, a  
78      INTEGER, SAVE:: nlev      INTEGER, SAVE:: nlev
79    
80      ! TEST SUR QSAT      ! TEST SUR QSAT
81      REAL p(iim + 1, jjm + 1, llmp1)      REAL p(iim + 1, jjm + 1, llmp1)
82      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)  
   
83      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.  
84    
85      !-----------------------------------------------------------------------      !-----------------------------------------------------------------------
86    
87      PRINT *, 'Call sequence information: guide'      !!PRINT *, 'Call sequence information: guide'
   
     ! calcul de l'humidite saturante  
88    
89      forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps      first_call: IF (itau == 0) THEN
90      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)  
91    
92      ! initialisations pour la lecture des reanalyses.         IF (online) THEN
93      ! 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  
94    
95      IF (online /= - 1) THEN            ! coordonnees du centre du zoom
96         IF (first) THEN            CALL coordij(clon, clat, ilon, ilat)
97            CALL conf_guide            ! aire de la maille au centre du zoom
98            file = 'guide'            aire_min = aire(ilon+(ilat - 1) * iip1)
99            CALL inigrads(igrads, rlonv, 180. / pi, -180., 180., rlatu, -90., 90., &            ! aire maximale de la maille
100                 180. / pi, presnivs, 1., dtgrads, file, 'dyn_zon ')            aire_max = 0.
101            PRINT *, '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'            DO ij = 1, ip1jmp1
102                 aire_max = max(aire_max, aire(ij))
103            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 ')  
104    
105               CALL wrgrads(igrads, llm, qsat, 'QSAT ', 'QSAT ')            factt = dtvr * iperiod / daysec
106    
107            END IF            CALL tau2alpha(3, factt, tau_min_v, tau_max_v, alpha_v)
108              CALL tau2alpha(2, factt, tau_min_u, tau_max_u, alpha_u)
109              CALL tau2alpha(1, factt, tau_min_t, tau_max_t, alpha_t)
110              CALL tau2alpha(1, factt, tau_min_q, tau_max_q, alpha_q)
111         ELSE         ELSE
112            count_no_rea = count_no_rea + 1            ! Cas où on force exactement par les variables analysées
113              alpha_t = 0.
114              alpha_u = 0.
115              alpha_v = 0.
116              alpha_q = 0.
117         END IF         END IF
118    
119         ! Guidage         step_rea = 1
120         ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses         count_no_rea = 0
121           ncid = -99
122         IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'  
123           ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
124         ditau = real(itau)         if (guide_u) call nf95_open('u.nc',Nf90_NOWRITe,ncid)
125         dday_step = real(day_step)         if (guide_v) call nf95_open('v.nc',nf90_nowrite,ncid)
126           if (guide_T) call nf95_open('T.nc',nf90_nowrite,ncid)
127         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  
128    
129         IF (guide_t) THEN         IF (ncep) THEN
130            DO l = 1, llm            call nf95_inq_dimid(ncid, 'LEVEL', dimid)
131               do j = 1, jjm + 1         ELSE
132                  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  
133         END IF         END IF
134           call nf95_inquire_dimension(ncid, dimid, nclen=nlev)
135         first = .FALSE.         PRINT *, 'nlev', nlev
136      end IF         call nf95_close(ncid)
137           ! Lecture du premier etat des reanalyses.
138           CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
139                masserea2, nlev)
140           qrea2 = max(qrea2, 0.1)
141    
142           CALL writefield("alpha_u", alpha_u)
143           CALL writefield("alpha_t", alpha_t)
144        END IF first_call
145    
146        ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
147    
148        ! Nudging fields are given 4 times per day:
149        IF (mod(itau, day_step / 4) == 0) THEN
150           vcovrea1 = vcovrea2
151           ucovrea1 = ucovrea2
152           tetarea1 = tetarea2
153           qrea1 = qrea2
154    
155           PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
156                count_no_rea, ' non lectures'
157           step_rea = step_rea + 1
158           CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
159                masserea2, nlev)
160           qrea2 = max(qrea2, 0.1)
161           factt = dtvr * iperiod / daysec
162           CALL writefield("ucov", ucov)
163           CALL writefield("ucovrea2", ucovrea2)
164           CALL writefield("teta", teta)
165           CALL writefield("tetarea2", tetarea2)
166           CALL writefield("qrea2", qrea2)
167           CALL writefield("q", q)
168        ELSE
169           count_no_rea = count_no_rea + 1
170        END IF
171    
172        ! Guidage
173    
174        tau = mod(real(itau) / real(day_step / 4), 1.)
175    
176        ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
177    
178        IF (guide_u) THEN
179           IF (itau == 0 .AND. ini_anal) then
180              ucov = ucovrea1
181           else
182              forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) * ucov(:, :, l) &
183                   + alpha_u * ((1. - tau) * ucovrea1(:, :, l) &
184                   + tau * ucovrea2(:, :, l))
185           end IF
186        END IF
187    
188        IF (guide_t) THEN
189           IF (itau == 0 .AND. ini_anal) then
190              teta = tetarea1
191           else
192              forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) * teta(:, :, l) &
193                   + alpha_t * ((1. - tau) * tetarea1(:, :, l) &
194                   + tau * tetarea2(:, :, l))
195           end IF
196        END IF
197    
198        IF (guide_q) THEN
199           ! Calcul de l'humidité saturante :
200           forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
201           CALL exner_hyb(ps, p, pks, pk)
202           qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
203    
204           ! humidité relative en % -> humidité spécifique
205           IF (itau == 0 .AND. ini_anal) then
206              q = qsat * qrea1 * 0.01
207           else
208              forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
209                   + alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
210                   + tau * qrea2(:, :, l)) * 0.01)
211           end IF
212        END IF
213    
214        IF (guide_v) THEN
215           IF (itau == 0 .AND. ini_anal) then
216              vcov = vcovrea1
217           else
218              forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) * vcov(:, :, l) &
219                   + alpha_v * ((1. - tau) * vcovrea1(:, :, l) &
220                   + tau * vcovrea2(:, :, l))
221           end IF
222        END IF
223    
224    END SUBROUTINE guide    END SUBROUTINE guide
225    
Legend:
Removed from v.91  
changed lines
  Added in v.109
  ViewVC Help
Powered by ViewVC 1.1.21