/[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/libf/dyn3d/guide.f90 revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/dyn3d/guide.f revision 112 by guez, Thu Sep 18 13:36:51 2014 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, version 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, 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  
   
   
   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  
   
     IMPLICIT NONE  
   
     !      ......   Version  du 10/01/98    ..........  
   
     !             avec  coordonnees  verticales hybrides  
     !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )  
   
     !=======================================================================  
     !  
     !   Auteur:  F.Hourdin  
     !   -------  
     !  
     !   Objet:  
     !   ------  
     !  
     !   GCM LMD nouvelle grille  
     !  
     !=======================================================================  
   
     !   ...  Dans inigeom , nouveaux calculs pour les elongations  cu , cv  
     !        et possibilite d'appeler une fonction f(y)  a derivee tangente  
     !        hyperbolique a la  place de la fonction a derivee sinusoidale.          
   
     !   ...  Possibilite de choisir le shema de Van-leer pour l'advection de  
     !         q  , en faisant iadv = 10  dans   traceur  (29/04/97) .  
     !  
     !-----------------------------------------------------------------------  
     !   Declarations:  
     !   -------------  
   
     include "netcdf.inc"  
7    
8      !   variables dynamiques    REAL aire_min, aire_max
     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./  
9    
10      print *,'Call sequence information: guide'  CONTAINS
11    
12      SUBROUTINE guide(itau, ucov, vcov, teta, q, ps)
13    
14        ! Author: F.Hourdin
15    
16        USE comconst, ONLY: cpp, daysec, dtvr, kappa
17        USE comgeom, ONLY: aire, rlatu, rlonv
18        USE conf_gcm_m, ONLY: day_step, iperiod
19        use conf_guide_m, only: conf_guide, guide_u, guide_v, guide_t, guide_q, &
20             ncep, ini_anal, tau_min_u, tau_max_u, tau_min_v, tau_max_v, &
21             tau_min_t, tau_max_t, tau_min_q, tau_max_q, online
22        USE dimens_m, ONLY: iim, jjm, llm
23        USE disvert_m, ONLY: ap, bp, preff, presnivs
24        USE exner_hyb_m, ONLY: exner_hyb
25        use netcdf, only: nf90_nowrite
26        use netcdf95, only: nf95_close, nf95_inq_dimid, nf95_inquire_dimension, &
27             nf95_open
28        use nr_util, only: pi
29        USE paramet_m, ONLY: iip1, ip1jmp1, jjp1, llmp1
30        USE q_sat_m, ONLY: q_sat
31        use read_reanalyse_m, only: read_reanalyse
32        USE serre, ONLY: clat, clon
33        use tau2alpha_m, only: tau2alpha
34        use writefield_m, only: writefield
35    
36        INTEGER, INTENT(IN):: itau
37        REAL, intent(inout):: ucov(:, :, :) ! (iim + 1, jjm + 1, llm) vent covariant
38        REAL, intent(inout):: vcov(:, :, :) ! (iim + 1, jjm, llm) ! vent covariant
39    
40        REAL, intent(inout):: teta(:, :, :) ! (iim + 1, jjm + 1, llm)
41        ! température potentielle
42    
43        REAL, intent(inout):: q(:, :, :) ! (iim + 1, jjm + 1, llm)
44        REAL, intent(in):: ps(:, :) ! (iim + 1, jjm + 1) pression au sol
45    
46        ! Local:
47    
48        ! variables dynamiques pour les réanalyses
49    
50        REAL, save:: ucovrea1(iim + 1, jjm + 1, llm), vcovrea1(iim + 1, jjm, llm)
51        ! vents covariants reanalyses
52    
53        REAL, save:: tetarea1(iim + 1, jjm + 1, llm) ! temp pot reales
54        REAL, save:: qrea1(iim + 1, jjm + 1, llm) ! temp pot reales
55    
56        REAL, save:: ucovrea2(iim + 1, jjm + 1, llm), vcovrea2(iim + 1, jjm, llm)
57        ! vents covariants reanalyses
58    
59        REAL, save:: tetarea2(iim + 1, jjm + 1, llm) ! temp pot reales
60        REAL, save:: qrea2(iim + 1, jjm + 1, llm) ! temp pot reales
61        REAL, save:: masserea2(ip1jmp1, llm) ! masse
62    
63        ! alpha determine la part des injections de donnees a chaque etape
64        ! alpha=0 signifie pas d'injection
65        ! alpha=1 signifie injection totale
66        REAL, save:: alpha_q(iim + 1, jjm + 1)
67        REAL, save:: alpha_t(iim + 1, jjm + 1)
68        REAL, save:: alpha_u(iim + 1, jjm + 1), alpha_v(iim + 1, jjm)
69    
70        INTEGER, save:: step_rea, count_no_rea
71    
72        INTEGER ilon, ilat
73        REAL factt ! pas de temps entre deux appels au guidage, en fraction de jour
74    
75        INTEGER ij, l
76        INTEGER ncid, dimid
77        REAL tau
78        INTEGER, SAVE:: nlev
79    
80        ! TEST SUR QSAT
81        REAL p(iim + 1, jjm + 1, llmp1)
82        real pk(iim + 1, jjm + 1, llm), pks(iim + 1, jjm + 1)
83        REAL qsat(iim + 1, jjm + 1, llm)
84    
85        !-----------------------------------------------------------------------
86    
87        !!PRINT *, 'Call sequence information: guide'
88    
89        first_call: IF (itau == 0) THEN
90           CALL conf_guide
91    
92           IF (online) THEN
93              ! Constantes de temps de rappel en jour
94    
95              ! coordonnees du centre du zoom
96              CALL coordij(clon, clat, ilon, ilat)
97              ! aire de la maille au centre du zoom
98              aire_min = aire(ilon+(ilat - 1) * iip1)
99              ! aire maximale de la maille
100              aire_max = 0.
101              DO ij = 1, ip1jmp1
102                 aire_max = max(aire_max, aire(ij))
103              END DO
104    
105              factt = dtvr * iperiod / daysec
106    
107              if (guide_u) CALL tau2alpha(3, factt, tau_min_v, tau_max_v, alpha_v)
108              if (guide_v) CALL tau2alpha(2, factt, tau_min_u, tau_max_u, alpha_u)
109              if (guide_t) CALL tau2alpha(1, factt, tau_min_t, tau_max_t, alpha_t)
110              if (guide_q) CALL tau2alpha(1, factt, tau_min_q, tau_max_q, alpha_q)
111           ELSE
112              ! Cas où on force exactement par les variables analysées
113              if (guide_u) alpha_t = 1.
114              if (guide_v) alpha_u = 1.
115              if (guide_t) alpha_v = 1.
116              if (guide_q) alpha_q = 1.
117           END IF
118    
119      !-----------------------------------------------------------------------         step_rea = 1
120      ! calcul de l'humidite saturante         count_no_rea = 0
     !-----------------------------------------------------------------------  
     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)  
121    
     !-----------------------------------------------------------------------  
   
     !-----------------------------------------------------------------------  
     !   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  
122         ! lecture d'un fichier netcdf pour determiner le nombre de niveaux         ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
123         if (guide_u) then         if (guide_u) then
124            if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)            call nf95_open('u.nc',Nf90_NOWRITe,ncid)
125         endif         else if (guide_v) then
126         !            call nf95_open('v.nc',nf90_nowrite,ncid)
127         if (guide_v) then         else if (guide_T) then
128            if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)            call nf95_open('T.nc',nf90_nowrite,ncid)
        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)  
129         else         else
130            status=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)            call nf95_open('hur.nc',nf90_nowrite, ncid)
131         endif         end if
        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:  
        !   ----------------------------------  
132    
133      endif ! first         IF (ncep) THEN
134      !            call nf95_inq_dimid(ncid, 'LEVEL', dimid)
135      !-----------------------------------------------------------------------         ELSE
136      !----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS:            call nf95_inq_dimid(ncid, 'PRESSURE', dimid)
137      !-----------------------------------------------------------------------         END IF
138           call nf95_inquire_dimension(ncid, dimid, nclen=nlev)
139           PRINT *, 'nlev', nlev
140           call nf95_close(ncid)
141           ! Lecture du premier etat des reanalyses.
142           CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
143                masserea2, nlev)
144           qrea2 = max(qrea2, 0.1)
145    
146           if (guide_u) CALL writefield("alpha_u", alpha_u)
147           if (guide_t) CALL writefield("alpha_t", alpha_t)
148        END IF first_call
149    
150        ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
151    
152        ! Nudging fields are given 4 times per day:
153        IF (mod(itau, day_step / 4) == 0) THEN
154           vcovrea1 = vcovrea2
155           ucovrea1 = ucovrea2
156           tetarea1 = tetarea2
157           qrea1 = qrea2
158    
159           PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
160                count_no_rea, ' non lectures'
161           step_rea = step_rea + 1
162           CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
163                masserea2, nlev)
164           qrea2 = max(qrea2, 0.1)
165           factt = dtvr * iperiod / daysec
166    
167      ditau=real(itau)         if (guide_u) then
168      DDAY_step=real(day_step)            CALL writefield("ucov", ucov)
169      write(*,*)'ditau,dday_step'            CALL writefield("ucovrea2", ucovrea2)
170      write(*,*)ditau,dday_step         end if
171      toto=4*ditau/dday_step  
172      reste=toto-aint(toto)         if (guide_t) then
173      !     write(*,*)'toto,reste',toto,reste            CALL writefield("teta", teta)
174              CALL writefield("tetarea2", tetarea2)
175      if (reste.eq.0.) then         end if
176         if (itau_test.eq.itau) then  
177            write(*,*)'deuxieme passage de advreel a itau=',itau         if (guide_q) then
178            stop            CALL writefield("qrea2", qrea2)
179              CALL writefield("q", q)
180           end if
181        ELSE
182           count_no_rea = count_no_rea + 1
183        END IF
184    
185        ! Guidage
186    
187        tau = mod(real(itau) / real(day_step / 4), 1.)
188    
189        ! x_gcm = a * x_gcm + (1 - a) * x_reanalyses
190    
191        IF (guide_u) THEN
192           IF (itau == 0 .AND. ini_anal) then
193              ucov = ucovrea1
194         else         else
195            vcovrea1(:,:)=vcovrea2(:,:)            forall (l = 1: llm) ucov(:, :, l) = (1. - alpha_u) * ucov(:, :, l) &
196            ucovrea1(:,:)=ucovrea2(:,:)                 + alpha_u * ((1. - tau) * ucovrea1(:, :, l) &
197            tetarea1(:,:)=tetarea2(:,:)                 + tau * ucovrea2(:, :, l))
198            qrea1(:,:)=qrea2(:,:)         end IF
199        END IF
200            print*,'LECTURE REANALYSES, pas ',step_rea &  
201                 ,'apres ',count_no_rea,' non lectures'      IF (guide_t) THEN
202            step_rea=step_rea+1         IF (itau == 0 .AND. ini_anal) then
203            itau_test=itau            teta = tetarea1
           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           ')  
   
     first=.false.  
   
     return  
   end subroutine guide  
   
   !=======================================================================  
   subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)  
     !=======================================================================  
   
     use dimens_m  
     use paramet_m  
     use comconst, only: pi  
     use comgeom  
     use serre  
     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))  
           enddo  
           zdx(1,j)=zdx(iip1,j)  
        enddo  
        do j=2,jjm  
           do i=1,iip1  
              zdy(i,j)=0.5*(cv_2d(i,j-1)+cv_2d(i,j))  
           enddo  
        enddo  
        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)  
        enddo  
        do j=1,jjp1  
           do i=1,iip1  
              dxdys(i,j)=sqrt(zdx(i,j)*zdx(i,j)+zdy(i,j)*zdy(i,j))  
           enddo  
        enddo  
        do j=1,jjp1  
           do i=1,iim  
              dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))  
           enddo  
           dxdyu(iip1,j)=dxdyu(1,j)  
        enddo  
        do j=1,jjm  
           do i=1,iip1  
              dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))  
           enddo  
        enddo  
   
        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))  
           enddo  
        enddo  
   
        if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then  
           print*,'ATTENTION modele peu zoome'  
           print*,'ATTENTION on prend une constante de guidage cste'  
           gamma=0.  
204         else         else
205            gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)            forall (l = 1: llm) teta(:, :, l) = (1. - alpha_t) * teta(:, :, l) &
206            print*,'gamma=',gamma                 + alpha_t * ((1. - tau) * tetarea1(:, :, l) &
207            if (gamma.lt.1.e-5) then                 + tau * tetarea2(:, :, l))
208               print*,'gamma =',gamma,'<1e-5'         end IF
209               stop      END IF
210            endif  
211            print*,'gamma=',gamma      IF (guide_q) THEN
212            gamma=log(0.5)/log(gamma)         ! Calcul de l'humidité saturante :
213         endif         forall (l = 1: llm + 1) p(:, :, l) = ap(l) + bp(l) * ps
214      endif         CALL exner_hyb(ps, p, pks, pk)
215           qsat = q_sat(pk * teta / cpp, preff * (pk / cpp)**(1. / kappa))
216      alphamin=factt/taumax  
217      alphamax=factt/taumin         ! humidité relative en % -> humidité spécifique
218           IF (itau == 0 .AND. ini_anal) then
219      do j=1,pjm            q = qsat * qrea1 * 0.01
220         do i=1,pim         else
221            if (type.eq.1) then            forall (l = 1: llm) q(:, :, l) = (1. - alpha_q) * q(:, :, l) &
222               dxdy_=dxdys(i,j)                 + alpha_q * (qsat(:, :, l) * ((1. - tau) * qrea1(:, :, l) &
223               zlat=rlatu(j)*180./pi                 + tau * qrea2(:, :, l)) * 0.01)
224            elseif (type.eq.2) then         end IF
225               dxdy_=dxdyu(i,j)      END IF
226               zlat=rlatu(j)*180./pi  
227            elseif (type.eq.3) then      IF (guide_v) THEN
228               dxdy_=dxdyv(i,j)         IF (itau == 0 .AND. ini_anal) then
229               zlat=rlatv(j)*180./pi            vcov = vcovrea1
230            endif         else
231            if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then            forall (l = 1: llm) vcov(:, :, l) = (1. - alpha_v) * vcov(:, :, l) &
232               !  pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin                 + alpha_v * ((1. - tau) * vcovrea1(:, :, l) &
233               alpha(i,j)=alphamin                 + tau * vcovrea2(:, :, l))
234            else         end IF
235               xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma      END IF
              xi=min(xi,1.)  
              if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then  
                 alpha(i,j)=xi*alphamin+(1.-xi)*alphamax  
              else  
                 alpha(i,j)=0.  
              endif  
           endif  
        enddo  
     enddo  
   
236    
237      return    END SUBROUTINE guide
   end subroutine tau2alpha  
238    
239  end module guide_m  END MODULE guide_m
Legend:
Removed from v.3  
changed lines
  Added in v.112
  ViewVC Help
Powered by ViewVC 1.1.21