/[lmdze]/trunk/Sources/dyn3d/Guide/guide.f
ViewVC logotype

Diff of /trunk/Sources/dyn3d/Guide/guide.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/dyn3d/guide.f revision 112 by guez, Thu Sep 18 13:36:51 2014 UTC trunk/Sources/dyn3d/Guide/guide.f revision 210 by guez, Tue Dec 13 16:02:23 2016 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, 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, &           tau_min_u, tau_max_u, tau_min_v, tau_max_v, tau_min_t, tau_max_t, &
18           ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &           tau_min_q, tau_max_q, online, factt
          tau_min_t, tau_max_t, tau_min_q, tau_max_q, online  
19      USE dimens_m, ONLY: iim, jjm, llm      USE dimens_m, ONLY: iim, jjm, llm
20      USE disvert_m, ONLY: ap, bp, preff, presnivs      USE disvert_m, ONLY: ap, bp, preff
21        use dynetat0_m, only: grossismx, grossismy, rlatu, rlatv
22      USE exner_hyb_m, ONLY: exner_hyb      USE exner_hyb_m, ONLY: exner_hyb
23      use netcdf, only: nf90_nowrite      use init_tau2alpha_m, only: init_tau2alpha
24      use netcdf95, only: nf95_close, nf95_inq_dimid, nf95_inquire_dimension, &      USE paramet_m, ONLY: iip1, jjp1
          nf95_open  
     use nr_util, only: pi  
     USE paramet_m, ONLY: iip1, ip1jmp1, jjp1, llmp1  
25      USE q_sat_m, ONLY: q_sat      USE q_sat_m, ONLY: q_sat
26      use read_reanalyse_m, only: read_reanalyse      use read_reanalyse_m, only: read_reanalyse
     USE serre, ONLY: clat, clon  
27      use tau2alpha_m, only: tau2alpha      use tau2alpha_m, only: tau2alpha
28      use writefield_m, only: writefield      use writefield_m, only: writefield
29    
# Line 45  CONTAINS Line 39  CONTAINS
39    
40      ! Local:      ! Local:
41    
42      ! variables dynamiques pour les réanalyses      ! Variables dynamiques pour les réanalyses
43    
44      REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)      REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
45      ! vents covariants reanalyses      ! vents covariants r\'eanalyses
46    
47      REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales      REAL, save:: tetarea1(iim + 1, jjm + 1, llm)
48      REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales      ! potential temperture from reanalysis
49        
50        REAL, save:: qrea1(iim + 1, jjm + 1, llm)
51    
52      REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)      REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
53      ! vents covariants reanalyses      ! vents covariants reanalyses
54    
55      REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales      REAL, save:: tetarea2(iim + 1, jjm + 1, llm)
56      REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales      ! potential temperture from reanalysis
57      REAL, save:: masserea2(ip1jmp1, llm) ! masse      
58        REAL, save:: qrea2(iim + 1, jjm + 1, llm)
59    
60      ! alpha determine la part des injections de donnees a chaque etape      ! alpha détermine la part des injections de données à chaque étape
61      ! alpha=0 signifie pas d'injection      ! alpha=0 signifie pas d'injection
62      ! alpha=1 signifie injection totale      ! alpha=1 signifie injection totale
63      REAL, save:: alpha_q(iim + 1, jjm + 1)      REAL, save:: alpha_q(iim + 1, jjm + 1)
64      REAL, save:: alpha_t(iim + 1, jjm + 1)      REAL, save:: alpha_t(iim + 1, jjm + 1)
65      REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)      REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)
66    
67      INTEGER, save:: step_rea, count_no_rea      INTEGER l
   
     INTEGER ilon, ilat  
     REAL factt ! pas de temps entre deux appels au guidage, en fraction de jour  
   
     INTEGER ij, l  
     INTEGER ncid, dimid  
68      REAL tau      REAL tau
     INTEGER, SAVE:: nlev  
69    
70      ! TEST SUR QSAT      ! TEST SUR QSAT
71      REAL p(iim + 1, jjm + 1, llmp1)      REAL p(iim + 1, jjm + 1, llm + 1)
72      real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)      real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
73      REAL qsat(iim + 1, jjm + 1, llm)      REAL qsat(iim + 1, jjm + 1, llm)
74    
75      !-----------------------------------------------------------------------      REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
   
     !!PRINT *, 'Call sequence information: guide'  
76    
77      first_call: IF (itau == 0) THEN      !-----------------------------------------------------------------------
        CALL conf_guide  
78    
79        IF (itau == 0) THEN
80         IF (online) THEN         IF (online) THEN
81            ! Constantes de temps de rappel en jour            IF (abs(grossismx - 1.) < 0.1 .OR. abs(grossismy - 1.) < 0.1) THEN
82                 ! grille regulière
83            ! coordonnees du centre du zoom               if (guide_u) alpha_u = 1. - exp(- factt / tau_max_u)
84            CALL coordij(clon, clat, ilon, ilat)               if (guide_v) alpha_v = 1. - exp(- factt / tau_max_v)
85            ! aire de la maille au centre du zoom               if (guide_t) alpha_t = 1. - exp(- factt / tau_max_t)
86            aire_min = aire(ilon+(ilat - 1) * iip1)               if (guide_q) alpha_q = 1. - exp(- factt / tau_max_q)
87            ! aire maximale de la maille            else
88            aire_max = 0.               call init_tau2alpha(dxdys, dxdyu, dxdyv)
89            DO ij = 1, ip1jmp1  
90               aire_max = max(aire_max, aire(ij))               if (guide_u) then
91            END DO                  CALL tau2alpha(dxdyu, rlatu, tau_min_u, tau_max_u, alpha_u)
92                    CALL writefield("alpha_u", alpha_u)
93            factt = dtvr * iperiod / daysec               end if
94    
95            if (guide_u) CALL tau2alpha(3, factt, tau_min_v, tau_max_v, alpha_v)               if (guide_v) then
96            if (guide_v) CALL tau2alpha(2, factt, tau_min_u, tau_max_u, alpha_u)                  CALL tau2alpha(dxdyv, rlatv, tau_min_v, tau_max_v, alpha_v)
97            if (guide_t) CALL tau2alpha(1, factt, tau_min_t, tau_max_t, alpha_t)                  CALL writefield("alpha_v", alpha_v)
98            if (guide_q) CALL tau2alpha(1, factt, tau_min_q, tau_max_q, alpha_q)               end if
99    
100                 if (guide_t) then
101                    CALL tau2alpha(dxdys, rlatu, tau_min_t, tau_max_t, alpha_t)
102                    CALL writefield("alpha_t", alpha_t)
103                 end if
104    
105                 if (guide_q)  then
106                    CALL tau2alpha(dxdys, rlatu, tau_min_q, tau_max_q, alpha_q)
107                    CALL writefield("alpha_q", alpha_q)
108                 end if
109              end IF
110         ELSE         ELSE
111            ! Cas où on force exactement par les variables analysées            ! Cas où on force exactement par les variables analysées
112            if (guide_u) alpha_t = 1.            if (guide_u) alpha_u = 1.
113            if (guide_v) alpha_u = 1.            if (guide_v) alpha_v = 1.
114            if (guide_t) alpha_v = 1.            if (guide_t) alpha_t = 1.
115            if (guide_q) alpha_q = 1.            if (guide_q) alpha_q = 1.
116         END IF         END IF
117    
118         step_rea = 1         ! Lecture du premier état des réanalyses :
119         count_no_rea = 0         CALL read_reanalyse(ps, ucovrea2, vcovrea2, tetarea2, qrea2)
   
        ! lecture d'un fichier netcdf pour determiner le nombre de niveaux  
        if (guide_u) then  
           call nf95_open('u.nc',Nf90_NOWRITe,ncid)  
        else if (guide_v) then  
           call nf95_open('v.nc',nf90_nowrite,ncid)  
        else if (guide_T) then  
           call nf95_open('T.nc',nf90_nowrite,ncid)  
        else  
           call nf95_open('hur.nc',nf90_nowrite, ncid)  
        end if  
   
        IF (ncep) THEN  
           call nf95_inq_dimid(ncid, 'LEVEL', dimid)  
        ELSE  
           call nf95_inq_dimid(ncid, 'PRESSURE', dimid)  
        END IF  
        call nf95_inquire_dimension(ncid, dimid, nclen=nlev)  
        PRINT *, 'nlev', nlev  
        call nf95_close(ncid)  
        ! Lecture du premier etat des reanalyses.  
        CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &  
             masserea2, nlev)  
120         qrea2 = max(qrea2, 0.1)         qrea2 = max(qrea2, 0.1)
121    
122         if (guide_u) CALL writefield("alpha_u", alpha_u)         if (ini_anal) then
123         if (guide_t) CALL writefield("alpha_t", alpha_t)            IF (guide_u) ucov = ucovrea2
124      END IF first_call            IF (guide_v) vcov = vcovrea2
125              IF (guide_t) teta = tetarea2
126    
127              IF (guide_q) then
128                 ! Calcul de l'humidité saturante :
129                 forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
130                 CALL exner_hyb(ps, p, pks, pk)
131                 q = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa)) &
132                      * qrea2 * 0.01
133              end IF
134           end if
135        END IF
136    
137      ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:      ! Importation des vents, pression et temp\'erature r\'eels :
138    
139      ! Nudging fields are given 4 times per day:      ! Nudging fields are given 4 times per day:
140      IF (mod(itau, day_step / 4) == 0) THEN      IF (mod(itau, day_step / 4) == 0) THEN
# Line 156  CONTAINS Line 143  CONTAINS
143         tetarea1 = tetarea2         tetarea1 = tetarea2
144         qrea1 = qrea2         qrea1 = qrea2
145    
146         PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &         CALL read_reanalyse(ps, ucovrea2, vcovrea2, tetarea2, qrea2)
             count_no_rea, ' non lectures'  
        step_rea = step_rea + 1  
        CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &  
             masserea2, nlev)  
147         qrea2 = max(qrea2, 0.1)         qrea2 = max(qrea2, 0.1)
        factt = dtvr * iperiod / daysec  
148    
149         if (guide_u) then         if (guide_u) then
150            CALL writefield("ucov", ucov)            CALL writefield("ucov", ucov)
# Line 178  CONTAINS Line 160  CONTAINS
160            CALL writefield("qrea2", qrea2)            CALL writefield("qrea2", qrea2)
161            CALL writefield("q", q)            CALL writefield("q", q)
162         end if         end if
     ELSE  
        count_no_rea = count_no_rea + 1  
163      END IF      END IF
164    
165      ! Guidage      ! Guidage
# Line 188  CONTAINS Line 168  CONTAINS
168    
169      ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses      ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
170    
171      IF (guide_u) THEN      IF (guide_u) forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) &
172         IF (itau == 0 .AND. ini_anal) then           * ucov(:, :, l) + alpha_u * ((1. - tau) * ucovrea1(:, :, l) + tau &
173            ucov = ucovrea1           * ucovrea2(:, :, l))
174         else  
175            forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) * ucov(:, :, l) &      IF (guide_v) forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) &
176                 + alpha_u * ((1. - tau) * ucovrea1(:, :, l) &           * vcov(:, :, l) + alpha_v * ((1. - tau) * vcovrea1(:, :, l) + tau &
177                 + tau * ucovrea2(:, :, l))           * vcovrea2(:, :, l))
178         end IF  
179      END IF      IF (guide_t) forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) &
180             * teta(:, :, l) + alpha_t * ((1. - tau) * tetarea1(:, :, l) + tau &
181      IF (guide_t) THEN           * tetarea2(:, :, l))
        IF (itau == 0 .AND. ini_anal) then  
           teta = tetarea1  
        else  
           forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) * teta(:, :, l) &  
                + alpha_t * ((1. - tau) * tetarea1(:, :, l) &  
                + tau * tetarea2(:, :, l))  
        end IF  
     END IF  
182    
183      IF (guide_q) THEN      IF (guide_q) THEN
184         ! Calcul de l'humidité saturante :         ! Calcul de l'humidité saturante :
# Line 215  CONTAINS Line 187  CONTAINS
187         qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))         qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
188    
189         ! humidité relative en % -> humidité spécifique         ! humidité relative en % -> humidité spécifique
190         IF (itau == 0 .AND. ini_anal) then         forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
191            q = qsat * qrea1 * 0.01              + alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
192         else              + tau * qrea2(:, :, l)) * 0.01)
           forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &  
                + alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &  
                + tau * qrea2(:, :, l)) * 0.01)  
        end IF  
     END IF  
   
     IF (guide_v) THEN  
        IF (itau == 0 .AND. ini_anal) then  
           vcov = vcovrea1  
        else  
           forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) * vcov(:, :, l) &  
                + alpha_v * ((1. - tau) * vcovrea1(:, :, l) &  
                + tau * vcovrea2(:, :, l))  
        end IF  
193      END IF      END IF
194    
195    END SUBROUTINE guide    END SUBROUTINE guide
Legend:
Removed from v.112  
changed lines
  Added in v.210
  ViewVC Help
Powered by ViewVC 1.1.21