/[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

trunk/dyn3d/guide.f revision 88 by guez, Tue Mar 11 15:09:02 2014 UTC trunk/dyn3d/Guide/guide.f revision 254 by guez, Mon Feb 5 10:39:38 2018 UTC
# Line 1  Line 1 
1  MODULE guide_m  MODULE guide_m
2    
3    ! From dyn3d/guide.F, version 1.3 2005/05/25 13:10:09    ! From dyn3d/guide.F, version 1.3, 2005/05/25 13:10:09
4    ! and dyn3d/guide.h, version 1.1.1.1 2004/05/19 12:53:06    ! and dyn3d/guide.h, version 1.1.1.1, 2004/05/19 12:53:06
5    
6    IMPLICIT NONE    IMPLICIT NONE
7    
   REAL aire_min, aire_max  
   
8  CONTAINS  CONTAINS
9    
10    SUBROUTINE guide(itau, ucov, vcov, teta, q, masse, ps)    SUBROUTINE guide(itau, ucov, vcov, teta, q, ps)
11    
12      ! Author: F.Hourdin      ! Author: F. Hourdin
13    
14      USE comconst, ONLY: cpp, daysec, dtvr, kappa      USE comconst, ONLY: cpp, kappa
15      USE comgeom, ONLY: aire, rlatu, rlonv      USE conf_gcm_m, ONLY: day_step
16      USE conf_gcm_m, ONLY: day_step, iperiod      use conf_guide_m, only: guide_u, guide_v, guide_t, guide_q, ini_anal, &
17      use conf_guide_m, only: conf_guide, guide_u, guide_v, guide_t, guide_q, &           alpha_u, alpha_v, alpha_t, alpha_q
          ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &  
          tau_min_t, tau_max_t, tau_min_q, tau_max_q, tau_min_p, tau_max_p, &  
          online  
18      USE dimens_m, ONLY: iim, jjm, llm      USE dimens_m, ONLY: iim, jjm, llm
19      USE disvert_m, ONLY: ap, bp, preff, presnivs      USE disvert_m, ONLY: ap, bp, preff
20      USE exner_hyb_m, ONLY: exner_hyb      USE exner_hyb_m, ONLY: exner_hyb
     USE inigrads_m, ONLY: inigrads  
     use massdair_m, only: massdair  
     use netcdf, only: nf90_nowrite, nf90_open, nf90_close, nf90_inq_dimid, &  
          nf90_inquire_dimension  
     use nr_util, only: pi  
     USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1, llmp1  
21      USE q_sat_m, ONLY: q_sat      USE q_sat_m, ONLY: q_sat
22      use read_reanalyse_m, only: read_reanalyse      use read_reanalyse_m, only: read_reanalyse
23      USE serre, ONLY: clat, clon      use writefield_m, only: writefield
     use tau2alpha_m, only: tau2alpha, dxdys  
24    
25      INTEGER, INTENT(IN):: itau      INTEGER, INTENT(IN):: itau
26        REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
27        REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
28    
29        REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
30        ! température potentielle
31    
32      ! variables dynamiques      REAL, intent(inout):: q(:, :, :) ! (iim + 1, jjm + 1, llm)
     REAL ucov(ip1jmp1, llm), vcov(ip1jm, llm) ! vents covariants  
     REAL, intent(inout):: teta(ip1jmp1, llm) ! temperature potentielle  
     REAL q(ip1jmp1, llm) ! temperature potentielle  
     REAL, intent(out):: masse(ip1jmp1, llm) ! masse d'air  
33      REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol      REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol
34    
35      ! Local:      ! Local:
36    
37      ! variables dynamiques pour les reanalyses.      ! Variables dynamiques pour les réanalyses
     REAL, save:: ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas  
     REAL, save:: tetarea1(ip1jmp1, llm) ! temp pot reales  
     REAL, save:: qrea1(ip1jmp1, llm) ! temp pot reales  
     REAL, save:: ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas  
     REAL, save:: tetarea2(ip1jmp1, llm) ! temp pot reales  
     REAL, save:: qrea2(ip1jmp1, llm) ! temp pot reales  
     REAL, save:: masserea2(ip1jmp1, llm) ! masse  
   
     REAL, save:: alpha_q(ip1jmp1)  
     REAL, save:: alpha_t(ip1jmp1), alpha_p(ip1jmp1)  
     REAL, save:: alpha_u(ip1jmp1), alpha_v(ip1jm)  
     REAL dday_step, toto, reste  
     real, save:: itau_test  
     INTEGER, save:: step_rea, count_no_rea  
   
     INTEGER ilon, ilat  
     REAL factt, ztau(ip1jmp1)  
   
     INTEGER ij, l  
     INTEGER ncidpl, varidpl, status  
     INTEGER rcod, rid  
     REAL ditau, tau, a  
     INTEGER, SAVE:: nlev  
38    
39      ! TEST SUR QSAT      REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
40      REAL p(iim + 1, jjm + 1, llmp1), pk(ip1jmp1, llm), pks(ip1jmp1)      ! vents covariants r\'eanalyses
     REAL pkf(ip1jmp1, llm)  
     REAL pres(ip1jmp1, llm)  
   
     REAL qsat(ip1jmp1, llm)  
     REAL unskap  
     REAL tnat(ip1jmp1, llm)  
   
     LOGICAL:: first = .TRUE.  
     CHARACTER(len=10) file  
     INTEGER:: igrads = 2  
     REAL:: dtgrads = 100.  
41    
42      !-----------------------------------------------------------------------      REAL, save:: tetarea1(iim + 1, jjm + 1, llm)
43        ! potential temperture from reanalysis
44    
45      PRINT *, 'Call sequence information: guide'      REAL, save:: qrea1(iim + 1, jjm + 1, llm)
46    
47      ! calcul de l'humidite saturante      REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
48        ! vents covariants reanalyses
49    
50      forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps      REAL, save:: tetarea2(iim + 1, jjm + 1, llm)
51      CALL massdair(p, masse)      ! potential temperture from reanalysis
     CALL exner_hyb(ps, p, pks, pk, pkf)  
     tnat(:, :) = pk(:, :)*teta(:, :)/cpp  
     unskap = 1./kappa  
     pres(:, :) = preff*(pk(:, :)/cpp)**unskap  
     qsat = q_sat(tnat, pres)  
   
     ! initialisations pour la lecture des reanalyses.  
     ! alpha determine la part des injections de donnees a chaque etape  
     ! alpha=1 signifie pas d'injection  
     ! alpha=0 signifie injection totale  
52    
53      IF (online==-1) THEN      REAL, save:: qrea2(iim + 1, jjm + 1, llm)
        RETURN  
     END IF  
54    
55      IF (first) THEN      INTEGER l
56         CALL conf_guide      REAL tau
        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) RETURN  
   
        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  
57    
58         ! itau_test montre si l'importation a deja ete faite au rang itau      ! TEST SUR QSAT
59         ! lecture d'un fichier netcdf pour determiner le nombre de niveaux      REAL p(iim + 1, jjm + 1, llm + 1)
60         if (guide_u) then      real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
61            if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)      REAL qsat(iim + 1, jjm + 1, llm)
62         endif  
63        !-----------------------------------------------------------------------
64    
65         if (guide_v) then      IF (itau == 0) THEN
66            if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)         ! Lecture du premier état des réanalyses :
67         endif         CALL read_reanalyse(ps, ucovrea2, vcovrea2, tetarea2, qrea2)
68           qrea2 = max(qrea2, 0.1)
69         if (guide_T) then  
70            if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)         if (ini_anal) then
71         endif            IF (guide_u) ucov = ucovrea2
72              IF (guide_v) vcov = vcovrea2
73         if (guide_Q) then            IF (guide_t) teta = tetarea2
74            if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)  
75         endif            IF (guide_q) then
76                 ! Calcul de l'humidité saturante :
77         IF (ncep) THEN               forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
78            status = nf90_inq_dimid(ncidpl, 'LEVEL', rid)               CALL exner_hyb(ps, p, pks, pk)
79         ELSE               q = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa)) &
80            status = nf90_inq_dimid(ncidpl, 'PRESSURE', rid)                    * qrea2 * 0.01
81         END IF            end IF
82         status = nf90_inquire_dimension(ncidpl, rid, len=nlev)         end if
83         PRINT *, 'nlev', nlev      END IF
84         rcod = nf90_close(ncidpl)  
85         ! Lecture du premier etat des reanalyses.      ! Importation des vents, pression et temp\'erature r\'eels :
86         CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &  
87              masserea2, nlev)      ! Nudging fields are given 4 times per day:
88         qrea2(:, :) = max(qrea2(:, :), 0.1)      IF (mod(itau, day_step / 4) == 0) THEN
89           vcovrea1 = vcovrea2
90         ! Debut de l'integration temporelle:         ucovrea1 = ucovrea2
91      END IF ! first         tetarea1 = tetarea2
92           qrea1 = qrea2
93      ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:  
94           CALL read_reanalyse(ps, ucovrea2, vcovrea2, tetarea2, qrea2)
95      ditau = real(itau)         qrea2 = max(qrea2, 0.1)
96      dday_step = real(day_step)  
97      WRITE (*, *) 'ditau, dday_step'         if (guide_u) then
98      WRITE (*, *) ditau, dday_step            CALL writefield("ucov", ucov)
99      toto = 4*ditau/dday_step            CALL writefield("ucovrea2", ucovrea2)
100      reste = toto - aint(toto)         end if
101    
102      IF (reste==0.) THEN         if (guide_t) then
103         IF (itau_test==itau) THEN            CALL writefield("teta", teta)
104            WRITE (*, *) 'deuxieme passage de advreel a itau=', itau            CALL writefield("tetarea2", tetarea2)
105            STOP         end if
106         ELSE  
107            vcovrea1(:, :) = vcovrea2(:, :)         if (guide_q) then
108            ucovrea1(:, :) = ucovrea2(:, :)            CALL writefield("qrea2", qrea2)
109            tetarea1(:, :) = tetarea2(:, :)            CALL writefield("q", q)
110            qrea1(:, :) = qrea2(:, :)         end if
   
           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(:), 1.E-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 ')  
   
           CALL wrgrads(igrads, llm, qsat, 'QSAT ', 'QSAT ')  
   
        END IF  
     ELSE  
        count_no_rea = count_no_rea + 1  
111      END IF      END IF
112    
113      ! Guidage      ! Guidage
     ! x_gcm = a * x_gcm + (1-a) * x_reanalyses  
114    
115      IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'      tau = mod(real(itau) / real(day_step / 4), 1.)
116    
117      ditau = real(itau)      ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
     dday_step = real(day_step)  
118    
119      tau = 4*ditau/dday_step      IF (guide_u) forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) &
120      tau = tau - aint(tau)           * ucov(:, :, l) + alpha_u * ((1. - tau) * ucovrea1(:, :, l) + tau &
121             * ucovrea2(:, :, l))
     ! 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  
122    
123      IF (guide_t) THEN      IF (guide_v) forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) &
124         DO l = 1, llm           * vcov(:, :, l) + alpha_v * ((1. - tau) * vcovrea1(:, :, l) + tau &
125            DO ij = 1, ip1jmp1           * vcovrea2(:, :, l))
              a = (1.-tau)*tetarea1(ij, l) + tau*tetarea2(ij, l)  
              teta(ij, l) = (1.-alpha_t(ij))*teta(ij, l) + alpha_t(ij)*a  
              IF (first .AND. ini_anal) teta(ij, l) = a  
           END DO  
        END DO  
     END IF  
126    
127      IF (guide_q) THEN      IF (guide_t) forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) &
128         DO l = 1, llm           * teta(:, :, l) + alpha_t * ((1. - tau) * tetarea1(:, :, l) + tau &
129            DO ij = 1, ip1jmp1           * tetarea2(:, :, l))
              a = (1.-tau)*qrea1(ij, l) + tau*qrea2(ij, l)  
              ! hum relative en % -> hum specif  
              a = qsat(ij, l)*a*0.01  
              q(ij, l) = (1.-alpha_q(ij))*q(ij, l) + alpha_q(ij)*a  
              IF (first .AND. ini_anal) q(ij, l) = a  
           END DO  
        END DO  
     END IF  
130    
131      ! vcov      IF (guide_q) THEN
132      IF (guide_v) THEN         ! Calcul de l'humidité saturante :
133         DO l = 1, llm         forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
134            DO ij = 1, ip1jm         CALL exner_hyb(ps, p, pks, pk)
135               a = (1.-tau)*vcovrea1(ij, l) + tau*vcovrea2(ij, l)         qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
136               vcov(ij, l) = (1.-alpha_v(ij))*vcov(ij, l) + alpha_v(ij)*a  
137               IF (first .AND. ini_anal) vcov(ij, l) = a         ! humidité relative en % -> humidité spécifique
138            END DO         forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
139            IF (first .AND. ini_anal) vcov(ij, l) = a              + alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
140         END DO              + tau * qrea2(:, :, l)) * 0.01)
141      END IF      END IF
142    
     first = .FALSE.  
   
143    END SUBROUTINE guide    END SUBROUTINE guide
144    
145  END MODULE guide_m  END MODULE guide_m

Legend:
Removed from v.88  
changed lines
  Added in v.254

  ViewVC Help
Powered by ViewVC 1.1.21