/[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/libf/dyn3d/guide.f90 revision 36 by guez, Thu Dec 2 17:11:04 2010 UTC trunk/dyn3d/guide.f revision 102 by guez, Tue Jul 15 13:43:24 2014 UTC
# Line 3  MODULE guide_m Line 3  MODULE guide_m
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    REAL tau_min_u, tau_max_u    IMPLICIT NONE
   REAL tau_min_v, tau_max_v  
   REAL tau_min_t, tau_max_t  
   REAL tau_min_q, tau_max_q  
   REAL tau_min_p, tau_max_p  
   REAL aire_min, aire_max  
   
7    
8    LOGICAL guide_u, guide_v, guide_t, guide_q, guide_p    REAL aire_min, aire_max
   REAL lat_min_guide, lat_max_guide  
   
   LOGICAL ncep, ini_anal  
   INTEGER online  
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    
16      USE dimens_m, ONLY : jjm, llm      USE comconst, ONLY: cpp, daysec, dtvr, kappa
17      USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1, llmp1      USE comgeom, ONLY: aire, rlatu, rlonv
18      USE comconst, ONLY : cpp, daysec, dtvr, kappa, pi      USE conf_gcm_m, ONLY: day_step, iperiod
19      USE comvert, ONLY : ap, bp, preff, presnivs      use conf_guide_m, only: conf_guide, guide_u, guide_v, guide_t, guide_q, &
20      USE conf_gcm_m, ONLY : day_step, iperiod           ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &
21      USE comgeom, ONLY : aire, rlatu, rlonv           tau_min_t, tau_max_t, tau_min_q, tau_max_q, tau_min_p, tau_max_p, &
22      USE serre, ONLY : clat, clon           online
23      USE q_sat_m, ONLY : q_sat      USE dimens_m, ONLY: iim, jjm, llm
24      USE exner_hyb_m, ONLY : exner_hyb      USE disvert_m, ONLY: ap, bp, preff, presnivs
25      USE pression_m, ONLY : pression      use dump2d_m, only: dump2d
26      USE inigrads_m, ONLY : inigrads      USE exner_hyb_m, ONLY: exner_hyb
27      use netcdf, only: nf90_nowrite, nf90_open, nf90_close      USE inigrads_m, ONLY: inigrads
28        use massdair_m, only: massdair
29      IMPLICIT NONE      use netcdf, only: nf90_nowrite, nf90_close, nf90_inq_dimid
30        use netcdf95, only: nf95_inquire_dimension, nf95_open
31      INCLUDE 'netcdf.inc'      use nr_util, only: pi
32        USE paramet_m, ONLY: iip1, ip1jm, ip1jmp1, jjp1, llmp1
33      !   variables dynamiques      USE q_sat_m, ONLY: q_sat
34      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants      use read_reanalyse_m, only: read_reanalyse
35      REAL, intent(inout):: teta(ip1jmp1, llm) ! temperature potentielle      USE serre, ONLY: clat, clon
36      REAL q(ip1jmp1, llm) ! temperature potentielle      use tau2alpha_m, only: tau2alpha, dxdys
37      REAL ps(ip1jmp1) ! pression  au sol  
38      REAL masse(ip1jmp1, llm) ! masse d'air      INTEGER, INTENT(IN):: itau
39    
40      !   common passe pour des sorties      ! variables dynamiques
41      REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)  
42      COMMON /comdxdy/dxdys, dxdyu, dxdyv      REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
43        REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
44      !   variables dynamiques pour les reanalyses.  
45      REAL ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas      REAL, intent(inout):: teta(iim + 1, jjm + 1, llm) ! température potentielle
46      REAL tetarea1(ip1jmp1, llm) ! temp pot  reales      REAL, intent(inout):: q(iim + 1, jjm + 1, llm)
47      REAL qrea1(ip1jmp1, llm) ! temp pot  reales      REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol
48      REAL psrea1(ip1jmp1) ! ps  
49      REAL ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas      ! Local:
50      REAL tetarea2(ip1jmp1, llm) ! temp pot  reales  
51      REAL qrea2(ip1jmp1, llm) ! temp pot  reales      ! variables dynamiques pour les reanalyses.
52      REAL masserea2(ip1jmp1, llm) ! masse  
53      REAL psrea2(ip1jmp1) ! ps      REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
54        ! vents covariants reanalyses
55      REAL alpha_q(ip1jmp1)  
56      REAL alpha_t(ip1jmp1), alpha_p(ip1jmp1)      REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales
57      REAL alpha_u(ip1jmp1), alpha_v(ip1jm)      REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales
58      REAL dday_step, toto, reste, itau_test  
59      INTEGER step_rea, count_no_rea      REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
60        ! vents covariants reanalyses
61    
62        REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales
63        REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales
64        REAL, save:: masserea2(ip1jmp1, llm) ! masse
65    
66        ! alpha determine la part des injections de donnees a chaque etape
67        ! alpha=1 signifie pas d'injection
68        ! alpha=0 signifie injection totale
69        REAL, save:: alpha_q(iim + 1, jjm + 1)
70        REAL, save:: alpha_t(iim + 1, jjm + 1), alpha_p(ip1jmp1)
71        REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)
72    
73        INTEGER, save:: step_rea, count_no_rea
74    
75      INTEGER ilon, ilat      INTEGER ilon, ilat
76      REAL factt, ztau(ip1jmp1)      REAL factt ! pas de temps entre deux appels au guidage, en fraction de jour
77        real ztau(iim + 1, jjm + 1)
78    
     INTEGER, INTENT (IN) :: itau  
79      INTEGER ij, l      INTEGER ij, l
80      INTEGER ncidpl, varidpl, nlev, status      INTEGER ncidpl, status
81      INTEGER rcod, rid      INTEGER rcod, rid
82      REAL ditau, tau, a      REAL tau
83      SAVE nlev      INTEGER, SAVE:: nlev
84    
85      !  TEST SUR QSAT      ! TEST SUR QSAT
86      REAL p(ip1jmp1, llmp1), pk(ip1jmp1, llm), pks(ip1jmp1)      REAL p(iim + 1, jjm + 1, llmp1)
87      REAL pkf(ip1jmp1, llm)      real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
     REAL pres(ip1jmp1, llm)  
   
     REAL qsat(ip1jmp1, llm)  
     REAL unskap  
     REAL tnat(ip1jmp1, llm)  
   
   
     LOGICAL first  
     SAVE first  
     DATA first/ .TRUE./  
   
     SAVE ucovrea1, vcovrea1, tetarea1, psrea1, qrea1  
     SAVE ucovrea2, vcovrea2, tetarea2, masserea2, psrea2, qrea2  
   
     SAVE alpha_t, alpha_q, alpha_u, alpha_v, alpha_p, itau_test  
     SAVE step_rea, count_no_rea  
   
     CHARACTER (10) file  
     INTEGER igrads  
     REAL dtgrads  
     SAVE igrads, dtgrads  
     DATA igrads, dtgrads/2, 100./  
88    
89      !-----------------------------------------------------------------------      REAL qsat(iim + 1, jjm + 1, llm)
90    
91      PRINT *, 'Call sequence information: guide'      INTEGER, parameter:: igrads = 2
92        REAL:: dtgrads = 100.
93    
94      ! calcul de l'humidite saturante      !-----------------------------------------------------------------------
   
     CALL pression(ip1jmp1, ap, bp, ps, p)  
     CALL massdair(p, masse)  
     PRINT *, 'OK1'  
     CALL exner_hyb(ps, p, pks, pk, pkf)  
     PRINT *, 'OK2'  
     tnat(:, :) = pk(:, :)*teta(:, :)/cpp  
     PRINT *, 'OK3'  
     unskap = 1./kappa  
     pres(:, :) = preff*(pk(:, :)/cpp)**unskap  
     PRINT *, 'OK4'  
     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  
   
     PRINT *, 'ONLINE=', online  
     IF (online==-1) THEN  
        RETURN  
     END IF  
95    
96      IF (first) THEN      PRINT *, 'Call sequence information: guide'
97    
98         PRINT *, 'initialisation du guide '      first_call: IF (itau == 0) THEN
99         CALL conf_guide         CALL conf_guide
100         PRINT *, 'apres conf_guide'         CALL inigrads(igrads, rlonv, 180. / pi, -180., 180., rlatu, -90., &
101                90., 180. / pi, presnivs, 1., dtgrads, 'guide', 'dyn_zon ')
        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)'  
102    
103         IF (online==-1) RETURN         IF (online) THEN
104         IF (online==1) THEN            ! Constantes de temps de rappel en jour
105    
106            !  Constantes de temps de rappel en jour            ! coordonnees du centre du zoom
           !  0.1 c'est en gros 2h30.  
           !  1e10  est une constante infinie donc en gros pas de guidage  
   
           !   coordonnees du centre du zoom  
107            CALL coordij(clon, clat, ilon, ilat)            CALL coordij(clon, clat, ilon, ilat)
108            !   aire de la maille au centre du zoom            ! aire de la maille au centre du zoom
109            aire_min = aire(ilon+(ilat-1)*iip1)            aire_min = aire(ilon+(ilat - 1) * iip1)
110            !   aire maximale de la maille            ! aire maximale de la maille
111            aire_max = 0.            aire_max = 0.
112            DO ij = 1, ip1jmp1            DO ij = 1, ip1jmp1
113               aire_max = max(aire_max, aire(ij))               aire_max = max(aire_max, aire(ij))
114            END DO            END DO
115            !  factt = pas de temps en fraction de jour  
116            factt = dtvr*iperiod/daysec            factt = dtvr * iperiod / daysec
117    
118            CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)            CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
119            CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)            CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
# Line 171  CONTAINS Line 122  CONTAINS
122            CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)            CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)
123    
124            CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')            CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
125            CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U   ')            CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U ')
126            CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T   ')            CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T ')
   
           !   Cas ou on force exactement par les variables analysees  
127         ELSE         ELSE
128              ! Cas ou on force exactement par les variables analysees
129            alpha_t = 0.            alpha_t = 0.
130            alpha_u = 0.            alpha_u = 0.
131            alpha_v = 0.            alpha_v = 0.
132            alpha_p = 0.            alpha_p = 0.
           !           physic=.false.  
133         END IF         END IF
134    
        itau_test = 1001  
135         step_rea = 1         step_rea = 1
136         count_no_rea = 0         count_no_rea = 0
137         ncidpl = -99         ncidpl = -99
138    
        !    itau_test    montre si l'importation a deja ete faite au rang itau  
139         ! lecture d'un fichier netcdf pour determiner le nombre de niveaux         ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
140         if (guide_u) then         if (guide_u) call nf95_open('u.nc',Nf90_NOWRITe,ncidpl)
141            if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)         if (guide_v) call nf95_open('v.nc',nf90_nowrite,ncidpl)
142         endif         if (guide_T) call nf95_open('T.nc',nf90_nowrite,ncidpl)
143           if (guide_Q) call nf95_open('hur.nc',nf90_nowrite, ncidpl)
        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  
144    
145         IF (ncep) THEN         IF (ncep) THEN
146            status = nf_inq_dimid(ncidpl, 'LEVEL', rid)            status = nf90_inq_dimid(ncidpl, 'LEVEL', rid)
147         ELSE         ELSE
148            status = nf_inq_dimid(ncidpl, 'PRESSURE', rid)            status = nf90_inq_dimid(ncidpl, 'PRESSURE', rid)
149         END IF         END IF
150         status = nf_inq_dimlen(ncidpl, rid, nlev)         call nf95_inquire_dimension(ncidpl, rid, nclen=nlev)
151         PRINT *, 'nlev', nlev         PRINT *, 'nlev', nlev
152         rcod = nf90_close(ncidpl)         rcod = nf90_close(ncidpl)
153         !   Lecture du premier etat des reanalyses.         ! Lecture du premier etat des reanalyses.
154         CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &         CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
155              masserea2, psrea2, 1, nlev)              masserea2, nlev)
156         qrea2(:, :) = max(qrea2(:, :), 0.1)         qrea2 = max(qrea2, 0.1)
157        END IF first_call
   
        !   Debut de l'integration temporelle:  
     END IF ! first  
158    
159      ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:      ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
160    
161      ditau = real(itau)      ! Nudging fields are given 4 times per day:
162      dday_step = real(day_step)      IF (mod(itau, day_step / 4) == 0) THEN
163      WRITE (*, *) 'ditau, dday_step'         vcovrea1 = vcovrea2
164      WRITE (*, *) ditau, dday_step         ucovrea1 = ucovrea2
165      toto = 4*ditau/dday_step         tetarea1 = tetarea2
166      reste = toto - aint(toto)         qrea1 = qrea2
167    
168      IF (reste==0.) THEN         PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
169         IF (itau_test==itau) THEN              count_no_rea, ' non lectures'
170            WRITE (*, *) 'deuxieme passage de advreel a itau=', itau         step_rea = step_rea + 1
171            STOP         CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
172         ELSE              masserea2, nlev)
173            vcovrea1(:, :) = vcovrea2(:, :)         qrea2 = max(qrea2, 0.1)
174            ucovrea1(:, :) = ucovrea2(:, :)         factt = dtvr * iperiod / daysec
175            tetarea1(:, :) = tetarea2(:, :)         ztau = factt / max(alpha_t, 1E-10)
176            qrea1(:, :) = qrea2(:, :)         CALL wrgrads(igrads, 1, aire, 'aire ', 'aire ')
177           CALL wrgrads(igrads, 1, dxdys, 'dxdy ', 'dxdy ')
178            PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &         CALL wrgrads(igrads, 1, alpha_u, 'au ', 'au ')
179                 count_no_rea, ' non lectures'         CALL wrgrads(igrads, 1, alpha_t, 'at ', 'at ')
180            step_rea = step_rea + 1         CALL wrgrads(igrads, 1, ztau, 'taut ', 'taut ')
181            itau_test = itau         CALL wrgrads(igrads, llm, ucov, 'u ', 'u ')
182            CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &         CALL wrgrads(igrads, llm, ucovrea2, 'ua ', 'ua ')
183                 qrea2, masserea2, psrea2, 1, nlev)         CALL wrgrads(igrads, llm, teta, 'T ', 'T ')
184            qrea2(:, :) = max(qrea2(:, :), 0.1)         CALL wrgrads(igrads, llm, tetarea2, 'Ta ', 'Ta ')
185            factt = dtvr*iperiod/daysec         CALL wrgrads(igrads, llm, qrea2, 'Qa ', 'Qa ')
186            ztau(:) = factt/max(alpha_t(:), 1.E-10)         CALL wrgrads(igrads, llm, q, 'Q ', 'Q ')
           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  
187      ELSE      ELSE
188         count_no_rea = count_no_rea + 1         count_no_rea = count_no_rea + 1
189      END IF      END IF
190    
191      !   Guidage      ! Guidage
     !    x_gcm = a * x_gcm + (1-a) * x_reanalyses  
   
     IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'  
192    
193      ditau = real(itau)      tau = mod(real(itau) / real(day_step / 4), 1.)
     dday_step = real(day_step)  
194    
195        ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
196    
     tau = 4*ditau/dday_step  
     tau = tau - aint(tau)  
   
     !  ucov  
197      IF (guide_u) THEN      IF (guide_u) THEN
198         DO l = 1, llm         IF (itau == 0 .AND. ini_anal) then
199            DO ij = 1, ip1jmp1            ucov = ucovrea1
200               a = (1.-tau)*ucovrea1(ij, l) + tau*ucovrea2(ij, l)         else
201               ucov(ij, l) = (1.-alpha_u(ij))*ucov(ij, l) + alpha_u(ij)*a            forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) * ucov(:, :, l) &
202               IF (first .AND. ini_anal) ucov(ij, l) = a                 + alpha_u * ((1. - tau) * ucovrea1(:, :, l) &
203            END DO                 + tau * ucovrea2(:, :, l))
204         END DO         end IF
205      END IF      END IF
206    
207      IF (guide_t) THEN      IF (guide_t) THEN
208         DO l = 1, llm         IF (itau == 0 .AND. ini_anal) then
209            DO ij = 1, ip1jmp1            teta = tetarea1
210               a = (1.-tau)*tetarea1(ij, l) + tau*tetarea2(ij, l)         else
211               teta(ij, l) = (1.-alpha_t(ij))*teta(ij, l) + alpha_t(ij)*a            forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) * teta(:, :, l) &
212               IF (first .AND. ini_anal) teta(ij, l) = a                 + alpha_t * ((1. - tau) * tetarea1(:, :, l) &
213            END DO                 + tau * tetarea2(:, :, l))
214         END DO         end IF
     END IF  
   
     !  P  
     IF (guide_p) THEN  
        DO ij = 1, ip1jmp1  
           a = (1.-tau)*psrea1(ij) + tau*psrea2(ij)  
           ps(ij) = (1.-alpha_p(ij))*ps(ij) + alpha_p(ij)*a  
           IF (first .AND. ini_anal) ps(ij) = a  
        END DO  
        CALL pression(ip1jmp1, ap, bp, ps, p)  
        CALL massdair(p, masse)  
215      END IF      END IF
216    
   
     !  q  
217      IF (guide_q) THEN      IF (guide_q) THEN
218         DO l = 1, llm         ! Calcul de l'humidité saturante :
219            DO ij = 1, ip1jmp1         forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
220               a = (1.-tau)*qrea1(ij, l) + tau*qrea2(ij, l)         CALL exner_hyb(ps, p, pks, pk)
221               !   hum relative en % -> hum specif         qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
222               a = qsat(ij, l)*a*0.01  
223               q(ij, l) = (1.-alpha_q(ij))*q(ij, l) + alpha_q(ij)*a         ! humidité relative en % -> humidité spécifique
224               IF (first .AND. ini_anal) q(ij, l) = a         IF (itau == 0 .AND. ini_anal) then
225            END DO            q = qsat * qrea1 * 0.01
226         END DO         else
227              forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
228                   + alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
229                   + tau * qrea2(:, :, l)) * 0.01)
230           end IF
231      END IF      END IF
232    
     ! vcov  
233      IF (guide_v) THEN      IF (guide_v) THEN
234         DO l = 1, llm         IF (itau == 0 .AND. ini_anal) then
235            DO ij = 1, ip1jm            vcov = vcovrea1
236               a = (1.-tau)*vcovrea1(ij, l) + tau*vcovrea2(ij, l)         else
237               vcov(ij, l) = (1.-alpha_v(ij))*vcov(ij, l) + alpha_v(ij)*a            forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) * vcov(:, :, l) &
238               IF (first .AND. ini_anal) vcov(ij, l) = a                 + alpha_v * ((1. - tau) * vcovrea1(:, :, l) &
239            END DO                 + tau * vcovrea2(:, :, l))
240            IF (first .AND. ini_anal) vcov(ij, l) = a         end IF
        END DO  
241      END IF      END IF
242    
     first = .FALSE.  
   
243    END SUBROUTINE guide    END SUBROUTINE guide
244    
   !=======================================================================  
   SUBROUTINE tau2alpha(type, pim, pjm, factt, taumin, taumax, alpha)  
     !=======================================================================  
   
     USE dimens_m, ONLY : iim, jjm  
     USE paramet_m, ONLY : iip1, jjp1  
     USE comconst, ONLY : pi  
     USE comgeom, ONLY : cu_2d, cv_2d, rlatu, rlatv  
     USE serre, ONLY : clat, clon, grossismx, grossismy  
     IMPLICIT NONE  
   
     !   arguments :  
     INTEGER type  
     INTEGER pim, pjm  
     REAL factt, taumin, taumax  
     REAL dxdy_, alpha(pim, pjm)  
     REAL dxdy_min, dxdy_max  
   
     !  local :  
     REAL alphamin, alphamax, gamma, xi  
     SAVE gamma  
     INTEGER i, j, ilon, ilat  
   
     LOGICAL first  
     SAVE first  
     DATA first/ .TRUE./  
   
     REAL zdx(iip1, jjp1), zdy(iip1, jjp1)  
   
     REAL zlat  
     REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)  
     COMMON /comdxdy/dxdys, dxdyu, dxdyv  
   
     IF (first) THEN  
        DO j = 2, jjm  
           DO i = 2, iip1  
              zdx(i, j) = 0.5*(cu_2d(i-1, j)+cu_2d(i, j))/cos(rlatu(j))  
           END DO  
           zdx(1, j) = zdx(iip1, j)  
        END DO  
        DO j = 2, jjm  
           DO i = 1, iip1  
              zdy(i, j) = 0.5*(cv_2d(i, j-1)+cv_2d(i, j))  
           END DO  
        END DO  
        DO i = 1, iip1  
           zdx(i, 1) = zdx(i, 2)  
           zdx(i, jjp1) = zdx(i, jjm)  
           zdy(i, 1) = zdy(i, 2)  
           zdy(i, jjp1) = zdy(i, jjm)  
        END DO  
        DO j = 1, jjp1  
           DO i = 1, iip1  
              dxdys(i, j) = sqrt(zdx(i, j)*zdx(i, j)+zdy(i, j)*zdy(i, j))  
           END DO  
        END DO  
        DO j = 1, jjp1  
           DO i = 1, iim  
              dxdyu(i, j) = 0.5*(dxdys(i, j)+dxdys(i+1, j))  
           END DO  
           dxdyu(iip1, j) = dxdyu(1, j)  
        END DO  
        DO j = 1, jjm  
           DO i = 1, iip1  
              dxdyv(i, j) = 0.5*(dxdys(i, j)+dxdys(i+1, j))  
           END DO  
        END DO  
   
        CALL dump2d(iip1, jjp1, dxdys, 'DX2DY2 SCAL  ')  
        CALL dump2d(iip1, jjp1, dxdyu, 'DX2DY2 U     ')  
        CALL dump2d(iip1, jjp1, dxdyv, 'DX2DY2 v     ')  
   
        !   coordonnees du centre du zoom  
        CALL coordij(clon, clat, ilon, ilat)  
        !   aire de la maille au centre du zoom  
        dxdy_min = dxdys(ilon, ilat)  
        !   dxdy maximale de la maille  
        dxdy_max = 0.  
        DO j = 1, jjp1  
           DO i = 1, iip1  
              dxdy_max = max(dxdy_max, dxdys(i, j))  
           END DO  
        END DO  
   
        IF (abs(grossismx-1.)<0.1 .OR. abs(grossismy-1.)<0.1) THEN  
           PRINT *, 'ATTENTION modele peu zoome'  
           PRINT *, 'ATTENTION on prend une constante de guidage cste'  
           gamma = 0.  
        ELSE  
           gamma = (dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)  
           PRINT *, 'gamma=', gamma  
           IF (gamma<1.E-5) THEN  
              PRINT *, 'gamma =', gamma, '<1e-5'  
              STOP  
           END IF  
           PRINT *, 'gamma=', gamma  
           gamma = log(0.5)/log(gamma)  
        END IF  
     END IF  
   
     alphamin = factt/taumax  
     alphamax = factt/taumin  
   
     DO j = 1, pjm  
        DO i = 1, pim  
           IF (type==1) THEN  
              dxdy_ = dxdys(i, j)  
              zlat = rlatu(j)*180./pi  
           ELSE IF (type==2) THEN  
              dxdy_ = dxdyu(i, j)  
              zlat = rlatu(j)*180./pi  
           ELSE IF (type==3) THEN  
              dxdy_ = dxdyv(i, j)  
              zlat = rlatv(j)*180./pi  
           END IF  
           IF (abs(grossismx-1.)<0.1 .OR. abs(grossismy-1.)<0.1) THEN  
              !  pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin  
              alpha(i, j) = alphamin  
           ELSE  
              xi = ((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma  
              xi = min(xi, 1.)  
              IF (lat_min_guide<=zlat .AND. zlat<=lat_max_guide) THEN  
                 alpha(i, j) = xi*alphamin + (1.-xi)*alphamax  
              ELSE  
                 alpha(i, j) = 0.  
              END IF  
           END IF  
        END DO  
     END DO  
   
   
     RETURN  
   END SUBROUTINE tau2alpha  
   
245  END MODULE guide_m  END MODULE guide_m

Legend:
Removed from v.36  
changed lines
  Added in v.102

  ViewVC Help
Powered by ViewVC 1.1.21