/[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 90 by guez, Wed Mar 12 21:16:36 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(iim + 1, jjm + 1, llm) ! température potentielle  
     REAL q(iim + 1, jjm + 1, llm)  
     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(iim + 1, jjm + 1, llm) ! temp pot reales  
     REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales  
     REAL, save:: ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas  
     REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales  
     REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales  
     REAL, save:: masserea2(ip1jmp1, llm) ! masse  
   
     REAL, save:: alpha_q(iim + 1, jjm + 1)  
     REAL, save:: alpha_t(iim + 1, jjm + 1), 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(iim + 1, jjm + 1)  
   
     INTEGER ij, i, j, l  
     INTEGER ncidpl, status  
     INTEGER rcod, rid  
     REAL ditau, tau, a  
     INTEGER, SAVE:: nlev  
   
     ! TEST SUR QSAT  
     REAL p(iim + 1, jjm + 1, llmp1)  
     real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)  
     REAL pres(iim + 1, jjm + 1, llm)  
38    
39      REAL qsat(iim + 1, jjm + 1, llm)      REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
40      REAL unskap      ! vents covariants r\'eanalyses
     REAL tnat(iim + 1, jjm + 1, llm)  
41    
42      LOGICAL:: first = .TRUE.      REAL, save:: tetarea1(iim + 1, jjm + 1, llm)
43      CHARACTER(len=10) file      ! potential temperture from reanalysis
     INTEGER:: igrads = 2  
     REAL:: dtgrads = 100.  
44    
45      !-----------------------------------------------------------------------      REAL, save:: qrea1(iim + 1, jjm + 1, llm)
46    
47      PRINT *, 'Call sequence information: guide'      REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
48        ! vents covariants reanalyses
49    
50      ! calcul de l'humidite saturante      REAL, save:: tetarea2(iim + 1, jjm + 1, llm)
51        ! potential temperture from reanalysis
52    
53      forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps      REAL, save:: qrea2(iim + 1, jjm + 1, llm)
     CALL massdair(p, masse)  
     CALL exner_hyb(ps, p, pks, pk)  
     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  
54    
55      IF (online== - 1) THEN      INTEGER l
56         RETURN      REAL tau
     END IF  
57    
58      IF (first) THEN      ! TEST SUR QSAT
59         CALL conf_guide      REAL p(iim + 1, jjm + 1, llm + 1)
60         file = 'guide'      real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
61         CALL inigrads(igrads, rlonv, 180. / pi, -180., 180., rlatu, -90., 90., &      REAL qsat(iim + 1, jjm + 1, llm)
             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  
62    
63         ! 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  
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)
   
        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)  
68         qrea2 = max(qrea2, 0.1)         qrea2 = max(qrea2, 0.1)
69    
70         ! Debut de l'integration temporelle:         if (ini_anal) then
71      END IF ! first            IF (guide_u) ucov = ucovrea2
72              IF (guide_v) vcov = vcovrea2
73              IF (guide_t) teta = tetarea2
74    
75              IF (guide_q) then
76                 ! Calcul de l'humidité saturante :
77                 forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
78                 CALL exner_hyb(ps, p, pks, pk)
79                 q = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa)) &
80                      * qrea2 * 0.01
81              end IF
82           end if
83        END IF
84    
85        ! Importation des vents, pression et temp\'erature r\'eels :
86    
87        ! Nudging fields are given 4 times per day:
88        IF (mod(itau, day_step / 4) == 0) THEN
89           vcovrea1 = vcovrea2
90           ucovrea1 = ucovrea2
91           tetarea1 = tetarea2
92           qrea1 = qrea2
93    
94      ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:         CALL read_reanalyse(ps, ucovrea2, vcovrea2, tetarea2, qrea2)
95           qrea2 = max(qrea2, 0.1)
96    
97      ditau = real(itau)         if (guide_u) then
98      dday_step = real(day_step)            CALL writefield("ucov", ucov)
99      WRITE (*, *) 'ditau, dday_step'            CALL writefield("ucovrea2", ucovrea2)
100      WRITE (*, *) ditau, dday_step         end if
101      toto = 4 * ditau / dday_step  
102      reste = toto - aint(toto)         if (guide_t) then
103              CALL writefield("teta", teta)
104      IF (reste==0.) THEN            CALL writefield("tetarea2", tetarea2)
105         IF (itau_test==itau) THEN         end if
106            WRITE (*, *) 'deuxieme passage de advreel a itau=', itau  
107            STOP         if (guide_q) then
108         ELSE            CALL writefield("qrea2", qrea2)
109            vcovrea1 = vcovrea2            CALL writefield("q", q)
110            ucovrea1 = ucovrea2         end if
           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 ')  
   
           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  
   
     IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'  
   
     ditau = real(itau)  
     dday_step = real(day_step)  
114    
115      tau = 4 * ditau / dday_step      tau = mod(real(itau) / real(day_step / 4), 1.)
     tau = tau - aint(tau)  
116    
117      ! ucov      ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
     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  
118    
119      IF (guide_t) THEN      IF (guide_u) forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) &
120         DO l = 1, llm           * ucov(:, :, l) + alpha_u * ((1. - tau) * ucovrea1(:, :, l) + tau &
121            do j = 1, jjm + 1           * ucovrea2(:, :, l))
122               DO i = 1, iim + 1  
123                  a = (1. - tau) * tetarea1(i, j, l) + tau * tetarea2(i, j, l)      IF (guide_v) forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) &
124                  teta(i, j, l) = (1. - alpha_t(i, j)) * teta(i, j, l) &           * vcov(:, :, l) + alpha_v * ((1. - tau) * vcovrea1(:, :, l) + tau &
125                       + alpha_t(i, j) * a           * vcovrea2(:, :, l))
126                  IF (first .AND. ini_anal) teta(i, j, l) = a  
127               END DO      IF (guide_t) forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) &
128            end do           * teta(:, :, l) + alpha_t * ((1. - tau) * tetarea1(:, :, l) + tau &
129         END DO           * tetarea2(:, :, l))
     END IF  
130    
131      IF (guide_q) THEN      IF (guide_q) THEN
132         DO l = 1, llm         ! Calcul de l'humidité saturante :
133            do j = 1, jjm + 1         forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
134               DO i = 1, iim + 1         CALL exner_hyb(ps, p, pks, pk)
135                  a = (1. - tau) * qrea1(i, j, l) + tau * qrea2(i, j, l)         qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
136                  ! hum relative en % -> hum specif  
137                  a = qsat(i, j, l) * a * 0.01         ! humidité relative en % -> humidité spécifique
138                  q(i, j, l) = (1. - alpha_q(i, j)) * q(i, j, l) &         forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
139                       + alpha_q(i, j) * a              + alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
140                  IF (first .AND. ini_anal) q(i, j, l) = a              + tau * qrea2(:, :, l)) * 0.01)
              END DO  
           end do  
        END DO  
141      END IF      END IF
142    
     ! 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  
     END IF  
   
     first = .FALSE.  
   
143    END SUBROUTINE guide    END SUBROUTINE guide
144    
145  END MODULE guide_m  END MODULE guide_m

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

  ViewVC Help
Powered by ViewVC 1.1.21