/[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 3 by guez, Wed Feb 27 13:16:39 2008 UTC revision 20 by guez, Wed Oct 15 16:19:57 2008 UTC
# Line 1  Line 1 
1  module guide_m  MODULE guide_m
2    
3    ! From dyn3d/guide.F,v 1.3 2005/05/25 13:10:09    ! From dyn3d/guide.F, v 1.3 2005/05/25 13:10:09
4    ! and dyn3d/guide.h,v 1.1.1.1 2004/05/19 12:53:06    ! and dyn3d/guide.h, v 1.1.1.1 2004/05/19 12:53:06
5    
6    real tau_min_u,tau_max_u    REAL :: tau_min_u, tau_max_u
7    real tau_min_v,tau_max_v    REAL :: tau_min_v, tau_max_v
8    real tau_min_T,tau_max_T    REAL :: tau_min_t, tau_max_t
9    real tau_min_q,tau_max_q    REAL :: tau_min_q, tau_max_q
10    real tau_min_P,tau_max_P    REAL :: tau_min_p, tau_max_p
11    real aire_min,aire_max    REAL :: aire_min, aire_max
   
   
   logical guide_u,guide_v,guide_T,guide_Q,guide_P  
   real lat_min_guide,lat_max_guide  
   
   LOGICAL ncep,ini_anal  
   integer online  
   
 contains  
   
   subroutine guide(itau,ucov,vcov,teta,q,masse,ps)  
   
     use dimens_m  
     use paramet_m  
     use comconst  
     use comdissnew  
     use comvert  
     use conf_gcm_m  
     use logic  
     use comgeom  
     use serre  
     use temps  
     use tracstoke  
     use ener  
     use q_sat_m, only: q_sat  
     use exner_hyb_m, only: exner_hyb  
     use pression_m, only: pression  
     use inigrads_m, only: inigrads  
12    
     IMPLICIT NONE  
13    
14      !      ......   Version  du 10/01/98    ..........    LOGICAL :: guide_u, guide_v, guide_t, guide_q, guide_p
15      REAL :: lat_min_guide, lat_max_guide
16    
17      !             avec  coordonnees  verticales hybrides    LOGICAL :: ncep, ini_anal
18      !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )    INTEGER :: online
19    
20      !=======================================================================  CONTAINS
21      !  
22      !   Auteur:  F.Hourdin  SUBROUTINE guide(itau, ucov, vcov, teta, q, masse, ps)
23      !   -------  
24      !    USE dimens_m, ONLY : jjm, llm
25      !   Objet:    USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1, llmp1
26      !   ------    USE comconst, ONLY : cpp, daysec, dtvr, kappa, pi
27      !    USE comvert, ONLY : ap, bp, preff, presnivs
28      !   GCM LMD nouvelle grille    USE conf_gcm_m, ONLY : day_step, iperiod
29      !    USE comgeom, ONLY : aire, rlatu, rlonv
30      !=======================================================================    USE serre, ONLY : clat, clon
31      USE q_sat_m, ONLY : q_sat
32      USE exner_hyb_m, ONLY : exner_hyb
33      USE pression_m, ONLY : pression
34      USE inigrads_m, ONLY : inigrads
35      use netcdf, only: nf90_nowrite, nf90_open, nf90_close
36    
37      IMPLICIT NONE
38      INCLUDE 'netcdf.inc'
39    
40      !      ......   Version  du 10/01/98    ..........
41    
42      !             avec  coordonnees  verticales hybrides
43      !   avec nouveaux operat. dissipation * ( gradiv2, divgrad2, nxgraro2 )
44    
45      !   ...  Dans inigeom , nouveaux calculs pour les elongations  cu , cv    !=======================================================================
46      !        et possibilite d'appeler une fonction f(y)  a derivee tangente  
47      !        hyperbolique a la  place de la fonction a derivee sinusoidale.            !   Auteur:  F.Hourdin
48      !   -------
49      !   ...  Possibilite de choisir le shema de Van-leer pour l'advection de  
50      !         q  , en faisant iadv = 10  dans   traceur  (29/04/97) .    !   Objet:
51      !    !   ------
52      !-----------------------------------------------------------------------  
53      !   Declarations:    !   GCM LMD nouvelle grille
     !   -------------  
   
     include "netcdf.inc"  
   
     !   variables dynamiques  
     REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants  
     REAL teta(ip1jmp1,llm)                 ! temperature potentielle  
     REAL q(ip1jmp1,llm)                 ! temperature potentielle  
     REAL ps(ip1jmp1)                       ! pression  au sol  
     REAL masse(ip1jmp1,llm)                ! masse d'air  
   
     !   common passe pour des sorties  
     real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)  
     common/comdxdy/dxdys,dxdyu,dxdyv  
   
     !   variables dynamiques pour les reanalyses.  
     REAL ucovrea1(ip1jmp1,llm),vcovrea1(ip1jm,llm) !vts cov reas  
     REAL tetarea1(ip1jmp1,llm)             ! temp pot  reales  
     REAL qrea1(ip1jmp1,llm)             ! temp pot  reales  
     REAL psrea1(ip1jmp1)             ! ps  
     REAL ucovrea2(ip1jmp1,llm),vcovrea2(ip1jm,llm) !vts cov reas  
     REAL tetarea2(ip1jmp1,llm)             ! temp pot  reales  
     REAL qrea2(ip1jmp1,llm)             ! temp pot  reales  
     REAL masserea2(ip1jmp1,llm)             ! masse  
     REAL psrea2(ip1jmp1)             ! ps  
   
     real alpha_q(ip1jmp1)  
     real alpha_T(ip1jmp1),alpha_P(ip1jmp1)  
     real alpha_u(ip1jmp1),alpha_v(ip1jm)  
     real dday_step,toto,reste,itau_test  
     INTEGER step_rea,count_no_rea  
   
     !IM 180305   real aire_min,aire_max  
     integer ilon,ilat  
     real factt,ztau(ip1jmp1)  
   
     INTEGER, intent(in):: itau  
     integer ij, l  
     integer ncidpl,varidpl,nlev,status  
     integer rcod,rid  
     real ditau,tau,a  
     save nlev  
   
     !  TEST SUR QSAT  
     real p(ip1jmp1,llmp1),pk(ip1jmp1,llm),pks(ip1jmp1)  
     real pkf(ip1jmp1,llm)  
     real pres(ip1jmp1,llm)  
   
     real qsat(ip1jmp1,llm)  
     real unskap  
     real tnat(ip1jmp1,llm)  
     !cccccccccccccccc  
   
   
     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./  
   
     print *,'Call sequence information: guide'  
   
     !-----------------------------------------------------------------------  
     ! 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.eq.-1) then  
        return  
     endif  
   
     if (first) then  
   
        print*,'initialisation du guide '  
        call conf_guide  
        print*,'apres conf_guide'  
   
        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.eq.-1) return  
        if (online.eq.1) then  
   
           !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
           !  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  
           !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
           !   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))  
           enddo  
           !  factt = pas de temps en fraction de jour  
           factt=dtvr*iperiod/daysec  
   
           !     subroutine tau2alpha(type,im,jm,factt,taumin,taumax,alpha)  
           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   ')  
   
           !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
           !   Cas ou on force exactement par les variables analysees  
        else  
           alpha_T=0.  
           alpha_u=0.  
           alpha_v=0.  
           alpha_P=0.  
           !           physic=.false.  
        endif  
   
        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) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)  
        endif  
        !  
        if (guide_v) then  
           if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)  
        endif  
        !  
        if (guide_T) then  
           if (ncidpl.eq.-99) ncidpl=NCOPN('T.nc',NCNOWRIT,rcod)  
        endif  
        !  
        if (guide_Q) then  
           if (ncidpl.eq.-99) ncidpl=NCOPN('hur.nc',NCNOWRIT,rcod)  
        endif  
        !  
        if (ncep) then  
           status=NF_INQ_DIMID(ncidpl,'LEVEL',rid)  
        else  
           status=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)  
        endif  
        status=NF_INQ_DIMLEN(ncidpl,rid,nlev)  
        print *,'nlev', nlev  
        call ncclos(ncidpl,rcod)  
        !   Lecture du premier etat des reanalyses.  
        call read_reanalyse(1,ps &  
             ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)  
        qrea2(:,:)=max(qrea2(:,:),0.1)  
   
   
        !-----------------------------------------------------------------------  
        !   Debut de l'integration temporelle:  
        !   ----------------------------------  
   
     endif ! 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)  
     !     write(*,*)'toto,reste',toto,reste  
   
     if (reste.eq.0.) then  
        if (itau_test.eq.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,psrea2,1,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      ' )  
   
        endif  
     else  
        count_no_rea=count_no_rea+1  
     endif  
   
     !-----------------------------------------------------------------------  
     !   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)  
   
   
     tau=4*ditau/dday_step  
     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  
           enddo  
        enddo  
     endif  
   
     !  teta  
     if (guide_T) then  
        do l=1,llm  
           do ij=1,ip1jmp1  
              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  
           enddo  
        enddo  
     endif  
   
     !  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  
        enddo  
        CALL pression(ip1jmp1,ap,bp,ps,p)  
        CALL massdair(p,masse)  
     endif  
   
   
     !  q  
     if (guide_Q) then  
        do l=1,llm  
           do ij=1,ip1jmp1  
              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  
           enddo  
        enddo  
     endif  
   
     ! 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  
           enddo  
           if (first.and.ini_anal) vcov(ij,l)=a  
        enddo  
     endif  
   
     !     call dump2d(iip1,jjp1,tetarea1,'TETA REA 1     ')  
     !     call dump2d(iip1,jjp1,tetarea2,'TETA REA 2     ')  
     !     call dump2d(iip1,jjp1,teta,'TETA           ')  
54    
55      first=.false.    !=======================================================================
56    
57      !   ...  Dans inigeom , nouveaux calculs pour les elongations  cu , cv
58      !        et possibilite d'appeler une fonction f(y)  a derivee tangente
59      !        hyperbolique a la  place de la fonction a derivee sinusoidale.        
60    
61      !   ...  Possibilite de choisir le shema de Van-leer pour l'advection de
62      !         q  , en faisant iadv = 10  dans   traceur  (29/04/97) .
63    
64      !-----------------------------------------------------------------------
65      !   Declarations:
66      !   -------------
67    
68    
69      !   variables dynamiques
70      REAL :: vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
71      REAL :: teta(ip1jmp1, llm) ! temperature potentielle
72      REAL :: q(ip1jmp1, llm) ! temperature potentielle
73      REAL :: ps(ip1jmp1) ! pression  au sol
74      REAL :: masse(ip1jmp1, llm) ! masse d'air
75    
76      !   common passe pour des sorties
77      REAL :: dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
78      COMMON /comdxdy/dxdys, dxdyu, dxdyv
79    
80      !   variables dynamiques pour les reanalyses.
81      REAL :: ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas
82      REAL :: tetarea1(ip1jmp1, llm) ! temp pot  reales
83      REAL :: qrea1(ip1jmp1, llm) ! temp pot  reales
84      REAL :: psrea1(ip1jmp1) ! ps
85      REAL :: ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas
86      REAL :: tetarea2(ip1jmp1, llm) ! temp pot  reales
87      REAL :: qrea2(ip1jmp1, llm) ! temp pot  reales
88      REAL :: masserea2(ip1jmp1, llm) ! masse
89      REAL :: psrea2(ip1jmp1) ! ps
90    
91      REAL :: alpha_q(ip1jmp1)
92      REAL :: alpha_t(ip1jmp1), alpha_p(ip1jmp1)
93      REAL :: alpha_u(ip1jmp1), alpha_v(ip1jm)
94      REAL :: dday_step, toto, reste, itau_test
95      INTEGER :: step_rea, count_no_rea
96    
97      !IM 180305   real aire_min, aire_max
98      INTEGER :: ilon, ilat
99      REAL :: factt, ztau(ip1jmp1)
100    
101      INTEGER, INTENT (IN) :: itau
102      INTEGER :: ij, l
103      INTEGER :: ncidpl, varidpl, nlev, status
104      INTEGER :: rcod, rid
105      REAL :: ditau, tau, a
106      SAVE nlev
107    
108      !  TEST SUR QSAT
109      REAL :: p(ip1jmp1, llmp1), pk(ip1jmp1, llm), pks(ip1jmp1)
110      REAL :: pkf(ip1jmp1, llm)
111      REAL :: pres(ip1jmp1, llm)
112    
113      REAL :: qsat(ip1jmp1, llm)
114      REAL :: unskap
115      REAL :: tnat(ip1jmp1, llm)
116      !cccccccccccccccc
117    
118    
119      LOGICAL :: first
120      SAVE first
121      DATA first/ .TRUE./
122    
123      SAVE ucovrea1, vcovrea1, tetarea1, psrea1, qrea1
124      SAVE ucovrea2, vcovrea2, tetarea2, masserea2, psrea2, qrea2
125    
126      SAVE alpha_t, alpha_q, alpha_u, alpha_v, alpha_p, itau_test
127      SAVE step_rea, count_no_rea
128    
129      CHARACTER (10) :: file
130      INTEGER :: igrads
131      REAL :: dtgrads
132      SAVE igrads, dtgrads
133      DATA igrads, dtgrads/2, 100./
134    
135      PRINT *, 'Call sequence information: guide'
136    
137      !-----------------------------------------------------------------------
138      ! calcul de l'humidite saturante
139      !-----------------------------------------------------------------------
140      CALL pression(ip1jmp1, ap, bp, ps, p)
141      CALL massdair(p, masse)
142      PRINT *, 'OK1'
143      CALL exner_hyb(ps, p, pks, pk, pkf)
144      PRINT *, 'OK2'
145      tnat(:, :) = pk(:, :)*teta(:, :)/cpp
146      PRINT *, 'OK3'
147      unskap = 1./kappa
148      pres(:, :) = preff*(pk(:, :)/cpp)**unskap
149      PRINT *, 'OK4'
150      qsat = q_sat(tnat, pres)
151    
152      !-----------------------------------------------------------------------
153    
154      !-----------------------------------------------------------------------
155      !   initialisations pour la lecture des reanalyses.
156      !    alpha determine la part des injections de donnees a chaque etape
157      !    alpha=1 signifie pas d'injection
158      !    alpha=0 signifie injection totale
159      !-----------------------------------------------------------------------
160    
161      PRINT *, 'ONLINE=', online
162      IF (online==-1) THEN
163         RETURN
164      END IF
165    
166      IF (first) THEN
167    
168         PRINT *, 'initialisation du guide '
169         CALL conf_guide
170         PRINT *, 'apres conf_guide'
171    
172         file = 'guide'
173         CALL inigrads(igrads, rlonv, 180./pi, -180., 180., rlatu, -90., 90., &
174              180./pi, presnivs, 1., dtgrads, file, 'dyn_zon ')
175    
176         PRINT *, &
177              '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
178    
179         IF (online==-1) RETURN
180         IF (online==1) THEN
181    
182            !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
183            !  Constantes de temps de rappel en jour
184            !  0.1 c'est en gros 2h30.
185            !  1e10  est une constante infinie donc en gros pas de guidage
186            !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
187            !   coordonnees du centre du zoom
188            CALL coordij(clon, clat, ilon, ilat)
189            !   aire de la maille au centre du zoom
190            aire_min = aire(ilon+(ilat-1)*iip1)
191            !   aire maximale de la maille
192            aire_max = 0.
193            DO ij = 1, ip1jmp1
194               aire_max = max(aire_max, aire(ij))
195            END DO
196            !  factt = pas de temps en fraction de jour
197            factt = dtvr*iperiod/daysec
198    
199            !     subroutine tau2alpha(type, im, jm, factt, taumin, taumax, alpha)
200            CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
201            CALL tau2alpha(2, iip1, jjp1, factt, tau_min_u, tau_max_u, alpha_u)
202            CALL tau2alpha(1, iip1, jjp1, factt, tau_min_t, tau_max_t, alpha_t)
203            CALL tau2alpha(1, iip1, jjp1, factt, tau_min_p, tau_max_p, alpha_p)
204            CALL tau2alpha(1, iip1, jjp1, factt, tau_min_q, tau_max_q, alpha_q)
205    
206            CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
207            CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U   ')
208            CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T   ')
209    
210            !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
211            !   Cas ou on force exactement par les variables analysees
212         ELSE
213            alpha_t = 0.
214            alpha_u = 0.
215            alpha_v = 0.
216            alpha_p = 0.
217            !           physic=.false.
218         END IF
219    
220         itau_test = 1001
221         step_rea = 1
222         count_no_rea = 0
223         ncidpl = -99
224    
225         !    itau_test    montre si l'importation a deja ete faite au rang itau
226         ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
227         if (guide_u) then
228            if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
229         endif
230    
231         if (guide_v) then
232            if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
233         endif
234    
235         if (guide_T) then
236            if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
237         endif
238    
239         if (guide_Q) then
240            if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
241         endif
242    
243         IF (ncep) THEN
244            status = nf_inq_dimid(ncidpl, 'LEVEL', rid)
245         ELSE
246            status = nf_inq_dimid(ncidpl, 'PRESSURE', rid)
247         END IF
248         status = nf_inq_dimlen(ncidpl, rid, nlev)
249         PRINT *, 'nlev', nlev
250         rcod = nf90_close(ncidpl)
251         !   Lecture du premier etat des reanalyses.
252         CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, masserea2, &
253              psrea2, 1, nlev)
254         qrea2(:, :) = max(qrea2(:, :), 0.1)
255    
256    
257         !-----------------------------------------------------------------------
258         !   Debut de l'integration temporelle:
259         !   ----------------------------------
260    
261      END IF ! first
262    
263      !-----------------------------------------------------------------------
264      !----- IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
265      !-----------------------------------------------------------------------
266    
267      ditau = real(itau)
268      dday_step = real(day_step)
269      WRITE (*, *) 'ditau, dday_step'
270      WRITE (*, *) ditau, dday_step
271      toto = 4*ditau/dday_step
272      reste = toto - aint(toto)
273      !     write(*, *)'toto, reste', toto, reste
274    
275      IF (reste==0.) THEN
276         IF (itau_test==itau) THEN
277            WRITE (*, *) 'deuxieme passage de advreel a itau=', itau
278            STOP
279         ELSE
280            vcovrea1(:, :) = vcovrea2(:, :)
281            ucovrea1(:, :) = ucovrea2(:, :)
282            tetarea1(:, :) = tetarea2(:, :)
283            qrea1(:, :) = qrea2(:, :)
284    
285            PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
286                 count_no_rea, ' non lectures'
287            step_rea = step_rea + 1
288            itau_test = itau
289            CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
290                 masserea2, psrea2, 1, nlev)
291            qrea2(:, :) = max(qrea2(:, :), 0.1)
292            factt = dtvr*iperiod/daysec
293            ztau(:) = factt/max(alpha_t(:), 1.E-10)
294            CALL wrgrads(igrads, 1, aire, 'aire      ', 'aire      ')
295            CALL wrgrads(igrads, 1, dxdys, 'dxdy      ', 'dxdy      ')
296            CALL wrgrads(igrads, 1, alpha_u, 'au        ', 'au        ')
297            CALL wrgrads(igrads, 1, alpha_t, 'at        ', 'at        ')
298            CALL wrgrads(igrads, 1, ztau, 'taut      ', 'taut      ')
299            CALL wrgrads(igrads, llm, ucov, 'u         ', 'u         ')
300            CALL wrgrads(igrads, llm, ucovrea2, 'ua        ', 'ua        ')
301            CALL wrgrads(igrads, llm, teta, 'T         ', 'T         ')
302            CALL wrgrads(igrads, llm, tetarea2, 'Ta        ', 'Ta        ')
303            CALL wrgrads(igrads, llm, qrea2, 'Qa        ', 'Qa        ')
304            CALL wrgrads(igrads, llm, q, 'Q         ', 'Q         ')
305    
306            CALL wrgrads(igrads, llm, qsat, 'QSAT      ', 'QSAT      ')
307    
308         END IF
309      ELSE
310         count_no_rea = count_no_rea + 1
311      END IF
312    
313      !-----------------------------------------------------------------------
314      !   Guidage
315      !    x_gcm = a * x_gcm + (1-a) * x_reanalyses
316      !-----------------------------------------------------------------------
317    
318      IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'
319    
320      ditau = real(itau)
321      dday_step = real(day_step)
322    
323    
324      tau = 4*ditau/dday_step
325      tau = tau - aint(tau)
326    
327      !  ucov
328      IF (guide_u) THEN
329         DO l = 1, llm
330            DO ij = 1, ip1jmp1
331               a = (1.-tau)*ucovrea1(ij, l) + tau*ucovrea2(ij, l)
332               ucov(ij, l) = (1.-alpha_u(ij))*ucov(ij, l) + alpha_u(ij)*a
333               IF (first .AND. ini_anal) ucov(ij, l) = a
334            END DO
335         END DO
336      END IF
337    
338      !  teta
339      IF (guide_t) THEN
340         DO l = 1, llm
341            DO ij = 1, ip1jmp1
342               a = (1.-tau)*tetarea1(ij, l) + tau*tetarea2(ij, l)
343               teta(ij, l) = (1.-alpha_t(ij))*teta(ij, l) + alpha_t(ij)*a
344               IF (first .AND. ini_anal) teta(ij, l) = a
345            END DO
346         END DO
347      END IF
348    
349      !  P
350      IF (guide_p) THEN
351         DO ij = 1, ip1jmp1
352            a = (1.-tau)*psrea1(ij) + tau*psrea2(ij)
353            ps(ij) = (1.-alpha_p(ij))*ps(ij) + alpha_p(ij)*a
354            IF (first .AND. ini_anal) ps(ij) = a
355         END DO
356         CALL pression(ip1jmp1, ap, bp, ps, p)
357         CALL massdair(p, masse)
358      END IF
359    
360    
361      !  q
362      IF (guide_q) THEN
363         DO l = 1, llm
364            DO ij = 1, ip1jmp1
365               a = (1.-tau)*qrea1(ij, l) + tau*qrea2(ij, l)
366               !   hum relative en % -> hum specif
367               a = qsat(ij, l)*a*0.01
368               q(ij, l) = (1.-alpha_q(ij))*q(ij, l) + alpha_q(ij)*a
369               IF (first .AND. ini_anal) q(ij, l) = a
370            END DO
371         END DO
372      END IF
373    
374      ! vcov
375      IF (guide_v) THEN
376         DO l = 1, llm
377            DO ij = 1, ip1jm
378               a = (1.-tau)*vcovrea1(ij, l) + tau*vcovrea2(ij, l)
379               vcov(ij, l) = (1.-alpha_v(ij))*vcov(ij, l) + alpha_v(ij)*a
380               IF (first .AND. ini_anal) vcov(ij, l) = a
381            END DO
382            IF (first .AND. ini_anal) vcov(ij, l) = a
383         END DO
384      END IF
385    
386      !     call dump2d(iip1, jjp1, tetarea1, 'TETA REA 1     ')
387      !     call dump2d(iip1, jjp1, tetarea2, 'TETA REA 2     ')
388      !     call dump2d(iip1, jjp1, teta, 'TETA           ')
389    
390      return    first = .FALSE.
391    end subroutine guide  
392      RETURN
393    END SUBROUTINE guide
394    
395    !=======================================================================    !=======================================================================
396    subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)    SUBROUTINE tau2alpha(type, pim, pjm, factt, taumin, taumax, alpha)
397      !=======================================================================      !=======================================================================
398    
399      use dimens_m      USE dimens_m, ONLY : iim, jjm
400      use paramet_m      USE paramet_m, ONLY : iip1, jjp1
401      use comconst, only: pi      USE comconst, ONLY : pi
402      use comgeom      USE comgeom, ONLY : cu_2d, cv_2d, rlatu, rlatv
403      use serre      USE serre, ONLY : clat, clon, grossismx, grossismy
404      implicit none      IMPLICIT NONE
405    
406      !   arguments :      !   arguments :
407      integer type      INTEGER :: type
408      integer pim,pjm      INTEGER :: pim, pjm
409      real factt,taumin,taumax      REAL :: factt, taumin, taumax
410      real dxdy_,alpha(pim,pjm)      REAL :: dxdy_, alpha(pim, pjm)
411      real dxdy_min,dxdy_max      REAL :: dxdy_min, dxdy_max
412    
413      !  local :      !  local :
414      real alphamin,alphamax,gamma,xi      REAL :: alphamin, alphamax, gamma, xi
415      save gamma      SAVE gamma
416      integer i,j,ilon,ilat      INTEGER :: i, j, ilon, ilat
417    
418      logical first      LOGICAL :: first
419      save first      SAVE first
420      data first/.true./      DATA first/ .TRUE./
421    
422      real zdx(iip1,jjp1),zdy(iip1,jjp1)      REAL :: zdx(iip1, jjp1), zdy(iip1, jjp1)
423    
424      real zlat      REAL :: zlat
425      real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)      REAL :: dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
426      common/comdxdy/dxdys,dxdyu,dxdyv      COMMON /comdxdy/dxdys, dxdyu, dxdyv
427    
428      if (first) then      IF (first) THEN
429         do j=2,jjm         DO j = 2, jjm
430            do i=2,iip1            DO i = 2, iip1
431               zdx(i,j)=0.5*(cu_2d(i-1,j)+cu_2d(i,j))/cos(rlatu(j))               zdx(i, j) = 0.5*(cu_2d(i-1, j)+cu_2d(i, j))/cos(rlatu(j))
432            enddo            END DO
433            zdx(1,j)=zdx(iip1,j)            zdx(1, j) = zdx(iip1, j)
434         enddo         END DO
435         do j=2,jjm         DO j = 2, jjm
436            do i=1,iip1            DO i = 1, iip1
437               zdy(i,j)=0.5*(cv_2d(i,j-1)+cv_2d(i,j))               zdy(i, j) = 0.5*(cv_2d(i, j-1)+cv_2d(i, j))
438            enddo            END DO
439         enddo         END DO
440         do i=1,iip1         DO i = 1, iip1
441            zdx(i,1)=zdx(i,2)            zdx(i, 1) = zdx(i, 2)
442            zdx(i,jjp1)=zdx(i,jjm)            zdx(i, jjp1) = zdx(i, jjm)
443            zdy(i,1)=zdy(i,2)            zdy(i, 1) = zdy(i, 2)
444            zdy(i,jjp1)=zdy(i,jjm)            zdy(i, jjp1) = zdy(i, jjm)
445         enddo         END DO
446         do j=1,jjp1         DO j = 1, jjp1
447            do i=1,iip1            DO i = 1, iip1
448               dxdys(i,j)=sqrt(zdx(i,j)*zdx(i,j)+zdy(i,j)*zdy(i,j))               dxdys(i, j) = sqrt(zdx(i, j)*zdx(i, j)+zdy(i, j)*zdy(i, j))
449            enddo            END DO
450         enddo         END DO
451         do j=1,jjp1         DO j = 1, jjp1
452            do i=1,iim            DO i = 1, iim
453               dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))               dxdyu(i, j) = 0.5*(dxdys(i, j)+dxdys(i+1, j))
454            enddo            END DO
455            dxdyu(iip1,j)=dxdyu(1,j)            dxdyu(iip1, j) = dxdyu(1, j)
456         enddo         END DO
457         do j=1,jjm         DO j = 1, jjm
458            do i=1,iip1            DO i = 1, iip1
459               dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))               dxdyv(i, j) = 0.5*(dxdys(i, j)+dxdys(i+1, j))
460            enddo            END DO
461         enddo         END DO
462    
463         call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL  ')         CALL dump2d(iip1, jjp1, dxdys, 'DX2DY2 SCAL  ')
464         call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U     ')         CALL dump2d(iip1, jjp1, dxdyu, 'DX2DY2 U     ')
465         call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v     ')         CALL dump2d(iip1, jjp1, dxdyv, 'DX2DY2 v     ')
466    
467         !   coordonnees du centre du zoom         !   coordonnees du centre du zoom
468         call coordij(clon,clat,ilon,ilat)         CALL coordij(clon, clat, ilon, ilat)
469         !   aire de la maille au centre du zoom         !   aire de la maille au centre du zoom
470         dxdy_min=dxdys(ilon,ilat)         dxdy_min = dxdys(ilon, ilat)
471         !   dxdy maximale de la maille         !   dxdy maximale de la maille
472         dxdy_max=0.         dxdy_max = 0.
473         do j=1,jjp1         DO j = 1, jjp1
474            do i=1,iip1            DO i = 1, iip1
475               dxdy_max=max(dxdy_max,dxdys(i,j))               dxdy_max = max(dxdy_max, dxdys(i, j))
476            enddo            END DO
477         enddo         END DO
478    
479         if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then         IF (abs(grossismx-1.)<0.1 .OR. abs(grossismy-1.)<0.1) THEN
480            print*,'ATTENTION modele peu zoome'            PRINT *, 'ATTENTION modele peu zoome'
481            print*,'ATTENTION on prend une constante de guidage cste'            PRINT *, 'ATTENTION on prend une constante de guidage cste'
482            gamma=0.            gamma = 0.
483         else         ELSE
484            gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)            gamma = (dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
485            print*,'gamma=',gamma            PRINT *, 'gamma=', gamma
486            if (gamma.lt.1.e-5) then            IF (gamma<1.E-5) THEN
487               print*,'gamma =',gamma,'<1e-5'               PRINT *, 'gamma =', gamma, '<1e-5'
488               stop               STOP
489            endif            END IF
490            print*,'gamma=',gamma            PRINT *, 'gamma=', gamma
491            gamma=log(0.5)/log(gamma)            gamma = log(0.5)/log(gamma)
492         endif         END IF
493      endif      END IF
494    
495      alphamin=factt/taumax      alphamin = factt/taumax
496      alphamax=factt/taumin      alphamax = factt/taumin
497    
498      do j=1,pjm      DO j = 1, pjm
499         do i=1,pim         DO i = 1, pim
500            if (type.eq.1) then            IF (type==1) THEN
501               dxdy_=dxdys(i,j)               dxdy_ = dxdys(i, j)
502               zlat=rlatu(j)*180./pi               zlat = rlatu(j)*180./pi
503            elseif (type.eq.2) then            ELSE IF (type==2) THEN
504               dxdy_=dxdyu(i,j)               dxdy_ = dxdyu(i, j)
505               zlat=rlatu(j)*180./pi               zlat = rlatu(j)*180./pi
506            elseif (type.eq.3) then            ELSE IF (type==3) THEN
507               dxdy_=dxdyv(i,j)               dxdy_ = dxdyv(i, j)
508               zlat=rlatv(j)*180./pi               zlat = rlatv(j)*180./pi
509            endif            END IF
510            if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then            IF (abs(grossismx-1.)<0.1 .OR. abs(grossismy-1.)<0.1) THEN
511               !  pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin               !  pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
512               alpha(i,j)=alphamin               alpha(i, j) = alphamin
513            else            ELSE
514               xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma               xi = ((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
515               xi=min(xi,1.)               xi = min(xi, 1.)
516               if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then               IF (lat_min_guide<=zlat .AND. zlat<=lat_max_guide) THEN
517                  alpha(i,j)=xi*alphamin+(1.-xi)*alphamax                  alpha(i, j) = xi*alphamin + (1.-xi)*alphamax
518               else               ELSE
519                  alpha(i,j)=0.                  alpha(i, j) = 0.
520               endif               END IF
521            endif            END IF
522         enddo         END DO
523      enddo      END DO
524    
525    
526      return      RETURN
527    end subroutine tau2alpha    END SUBROUTINE tau2alpha
528    
529  end module guide_m  END MODULE guide_m

Legend:
Removed from v.3  
changed lines
  Added in v.20

  ViewVC Help
Powered by ViewVC 1.1.21