/[lmdze]/trunk/libf/dyn3d/guide.f90
ViewVC logotype

Diff of /trunk/libf/dyn3d/guide.f90

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

revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC revision 37 by guez, Tue Dec 21 15:45:48 2010 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:  
     !   -------------  
7    
8      include "netcdf.inc"    REAL tau_min_u, tau_max_u
9      REAL tau_min_v, tau_max_v
10      REAL tau_min_t, tau_max_t
11      REAL tau_min_q, tau_max_q
12      REAL tau_min_p, tau_max_p
13      REAL aire_min, aire_max
14    
     !   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  
15    
16      !   common passe pour des sorties    LOGICAL guide_u, guide_v, guide_t, guide_q, guide_p
17      real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)    LOGICAL ncep, ini_anal
18      common/comdxdy/dxdys,dxdyu,dxdyv    INTEGER online
19    
20      !   variables dynamiques pour les reanalyses.  CONTAINS
     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  
21    
22      !  TEST SUR QSAT    SUBROUTINE guide(itau, ucov, vcov, teta, q, masse, ps)
     real p(ip1jmp1,llmp1),pk(ip1jmp1,llm),pks(ip1jmp1)  
     real pkf(ip1jmp1,llm)  
     real pres(ip1jmp1,llm)  
23    
24      real qsat(ip1jmp1,llm)      ! Author: F.Hourdin
     real unskap  
     real tnat(ip1jmp1,llm)  
     !cccccccccccccccc  
25    
26        USE dimens_m, ONLY : jjm, llm
27        USE paramet_m, ONLY : iip1, ip1jm, ip1jmp1, jjp1, llmp1
28        USE comconst, ONLY : cpp, daysec, dtvr, kappa, pi
29        USE comvert, ONLY : ap, bp, preff, presnivs
30        USE conf_gcm_m, ONLY : day_step, iperiod
31        USE comgeom, ONLY : aire, rlatu, rlonv
32        USE serre, ONLY : clat, clon
33        USE q_sat_m, ONLY : q_sat
34        USE exner_hyb_m, ONLY : exner_hyb
35        USE inigrads_m, ONLY : inigrads
36        use netcdf, only: nf90_nowrite, nf90_open, nf90_close
37        use tau2alpha_m, only: tau2alpha
38    
39      LOGICAL first      INCLUDE 'netcdf.inc'
     save first  
     data first/.true./  
40    
41      save ucovrea1,vcovrea1,tetarea1,psrea1,qrea1      !   variables dynamiques
42      save ucovrea2,vcovrea2,tetarea2,masserea2,psrea2,qrea2      REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
43        REAL, intent(inout):: teta(ip1jmp1, llm) ! temperature potentielle
44        REAL q(ip1jmp1, llm) ! temperature potentielle
45        REAL ps(ip1jmp1) ! pression  au sol
46        REAL masse(ip1jmp1, llm) ! masse d'air
47    
48      save alpha_T,alpha_q,alpha_u,alpha_v,alpha_P,itau_test      !   common passe pour des sorties
49      save step_rea,count_no_rea      REAL dxdys(iip1, jjp1), dxdyu(iip1, jjp1), dxdyv(iip1, jjm)
50        COMMON /comdxdy/dxdys, dxdyu, dxdyv
51    
52      character*10 file      !   variables dynamiques pour les reanalyses.
53      integer igrads      REAL ucovrea1(ip1jmp1, llm), vcovrea1(ip1jm, llm) !vts cov reas
54      real dtgrads      REAL tetarea1(ip1jmp1, llm) ! temp pot  reales
55      save igrads,dtgrads      REAL qrea1(ip1jmp1, llm) ! temp pot  reales
56      data igrads,dtgrads/2,100./      REAL psrea1(ip1jmp1) ! ps
57        REAL ucovrea2(ip1jmp1, llm), vcovrea2(ip1jm, llm) !vts cov reas
58        REAL tetarea2(ip1jmp1, llm) ! temp pot  reales
59        REAL qrea2(ip1jmp1, llm) ! temp pot  reales
60        REAL masserea2(ip1jmp1, llm) ! masse
61        REAL psrea2(ip1jmp1) ! ps
62    
63        REAL alpha_q(ip1jmp1)
64        REAL alpha_t(ip1jmp1), alpha_p(ip1jmp1)
65        REAL alpha_u(ip1jmp1), alpha_v(ip1jm)
66        REAL dday_step, toto, reste, itau_test
67        INTEGER step_rea, count_no_rea
68    
69        INTEGER ilon, ilat
70        REAL factt, ztau(ip1jmp1)
71    
72        INTEGER, INTENT (IN) :: itau
73        INTEGER ij, l
74        INTEGER ncidpl, varidpl, nlev, status
75        INTEGER rcod, rid
76        REAL ditau, tau, a
77        SAVE nlev
78    
79      print *,'Call sequence information: guide'      !  TEST SUR QSAT
80        REAL p(ip1jmp1, llmp1), pk(ip1jmp1, llm), pks(ip1jmp1)
81        REAL pkf(ip1jmp1, llm)
82        REAL pres(ip1jmp1, llm)
83    
84        REAL qsat(ip1jmp1, llm)
85        REAL unskap
86        REAL tnat(ip1jmp1, llm)
87    
88        LOGICAL:: first = .TRUE.
89    
90        SAVE ucovrea1, vcovrea1, tetarea1, psrea1, qrea1
91        SAVE ucovrea2, vcovrea2, tetarea2, masserea2, psrea2, qrea2
92    
93        SAVE alpha_t, alpha_q, alpha_u, alpha_v, alpha_p, itau_test
94        SAVE step_rea, count_no_rea
95    
96        CHARACTER (10) file
97        INTEGER igrads
98        REAL dtgrads
99        SAVE igrads, dtgrads
100        DATA igrads, dtgrads/2, 100./
101    
102      !-----------------------------------------------------------------------      !-----------------------------------------------------------------------
103    
104        PRINT *, 'Call sequence information: guide'
105    
106      ! calcul de l'humidite saturante      ! 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)  
107    
108      !-----------------------------------------------------------------------      forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
109        CALL massdair(p, masse)
110        PRINT *, 'OK1'
111        CALL exner_hyb(ps, p, pks, pk, pkf)
112        PRINT *, 'OK2'
113        tnat(:, :) = pk(:, :)*teta(:, :)/cpp
114        PRINT *, 'OK3'
115        unskap = 1./kappa
116        pres(:, :) = preff*(pk(:, :)/cpp)**unskap
117        PRINT *, 'OK4'
118        qsat = q_sat(tnat, pres)
119    
     !-----------------------------------------------------------------------  
120      !   initialisations pour la lecture des reanalyses.      !   initialisations pour la lecture des reanalyses.
121      !    alpha determine la part des injections de donnees a chaque etape      !    alpha determine la part des injections de donnees a chaque etape
122      !    alpha=1 signifie pas d'injection      !    alpha=1 signifie pas d'injection
123      !    alpha=0 signifie injection totale      !    alpha=0 signifie injection totale
     !-----------------------------------------------------------------------  
124    
125      print*,'ONLINE=',online      PRINT *, 'ONLINE=', online
126      if(online.eq.-1) then      IF (online==-1) THEN
127         return         RETURN
128      endif      END IF
129    
130      if (first) then      IF (first) THEN
131    
132         print*,'initialisation du guide '         PRINT *, 'initialisation du guide '
133         call conf_guide         CALL conf_guide
134         print*,'apres conf_guide'         PRINT *, 'apres conf_guide'
135    
136         file='guide'         file = 'guide'
137         call inigrads(igrads &         CALL inigrads(igrads, rlonv, 180./pi, -180., 180., rlatu, -90., 90., &
138              ,rlonv,180./pi,-180.,180.,rlatu,-90.,90.,180./pi &              180./pi, presnivs, 1., dtgrads, file, 'dyn_zon ')
             ,presnivs,1. &  
             ,dtgrads,file,'dyn_zon ')  
139    
140         print* &         PRINT *, &
141              ,'1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'              '1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
142    
143         if(online.eq.-1) return         IF (online==-1) RETURN
144         if (online.eq.1) then         IF (online==1) THEN
145    
           !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
146            !  Constantes de temps de rappel en jour            !  Constantes de temps de rappel en jour
147            !  0.1 c'est en gros 2h30.            !  0.1 c'est en gros 2h30.
148            !  1e10  est une constante infinie donc en gros pas de guidage            !  1e10  est une constante infinie donc en gros pas de guidage
149            !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
150            !   coordonnees du centre du zoom            !   coordonnees du centre du zoom
151            call coordij(clon,clat,ilon,ilat)            CALL coordij(clon, clat, ilon, ilat)
152            !   aire de la maille au centre du zoom            !   aire de la maille au centre du zoom
153            aire_min=aire(ilon+(ilat-1)*iip1)            aire_min = aire(ilon+(ilat-1)*iip1)
154            !   aire maximale de la maille            !   aire maximale de la maille
155            aire_max=0.            aire_max = 0.
156            do ij=1,ip1jmp1            DO ij = 1, ip1jmp1
157               aire_max=max(aire_max,aire(ij))               aire_max = max(aire_max, aire(ij))
158            enddo            END DO
159            !  factt = pas de temps en fraction de jour            !  factt = pas de temps en fraction de jour
160            factt=dtvr*iperiod/daysec            factt = dtvr*iperiod/daysec
161    
162            !     subroutine tau2alpha(type,im,jm,factt,taumin,taumax,alpha)            CALL tau2alpha(3, iip1, jjm, factt, tau_min_v, tau_max_v, alpha_v)
163            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)
164            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)
165            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)
166            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)
167            call tau2alpha(1,iip1,jjp1,factt,tau_min_q,tau_max_q,alpha_q)  
168              CALL dump2d(iip1, jjp1, aire, 'AIRE MAILLe ')
169            call dump2d(iip1,jjp1,aire,'AIRE MAILLe ')            CALL dump2d(iip1, jjp1, alpha_u, 'COEFF U   ')
170            call dump2d(iip1,jjp1,alpha_u,'COEFF U   ')            CALL dump2d(iip1, jjp1, alpha_t, 'COEFF T   ')
           call dump2d(iip1,jjp1,alpha_T,'COEFF T   ')  
171    
           !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
172            !   Cas ou on force exactement par les variables analysees            !   Cas ou on force exactement par les variables analysees
173         else         ELSE
174            alpha_T=0.            alpha_t = 0.
175            alpha_u=0.            alpha_u = 0.
176            alpha_v=0.            alpha_v = 0.
177            alpha_P=0.            alpha_p = 0.
178            !           physic=.false.            !           physic=.false.
179         endif         END IF
180    
181         itau_test=1001         itau_test = 1001
182         step_rea=1         step_rea = 1
183         count_no_rea=0         count_no_rea = 0
184         ncidpl=-99         ncidpl = -99
185    
186         !    itau_test    montre si l'importation a deja ete faite au rang itau         !    itau_test    montre si l'importation a deja ete faite au rang itau
187         ! lecture d'un fichier netcdf pour determiner le nombre de niveaux         ! lecture d'un fichier netcdf pour determiner le nombre de niveaux
188         if (guide_u) then         if (guide_u) then
189            if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)            if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
190         endif         endif
191         !  
192         if (guide_v) then         if (guide_v) then
193            if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)            if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
194         endif         endif
195         !  
196         if (guide_T) then         if (guide_T) then
197            if (ncidpl.eq.-99) ncidpl=NCOPN('T.nc',NCNOWRIT,rcod)            if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
198         endif         endif
199         !  
200         if (guide_Q) then         if (guide_Q) then
201            if (ncidpl.eq.-99) ncidpl=NCOPN('hur.nc',NCNOWRIT,rcod)            if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
202         endif         endif
203         !  
204         if (ncep) then         IF (ncep) THEN
205            status=NF_INQ_DIMID(ncidpl,'LEVEL',rid)            status = nf_inq_dimid(ncidpl, 'LEVEL', rid)
206         else         ELSE
207            status=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)            status = nf_inq_dimid(ncidpl, 'PRESSURE', rid)
208         endif         END IF
209         status=NF_INQ_DIMLEN(ncidpl,rid,nlev)         status = nf_inq_dimlen(ncidpl, rid, nlev)
210         print *,'nlev', nlev         PRINT *, 'nlev', nlev
211         call ncclos(ncidpl,rcod)         rcod = nf90_close(ncidpl)
212         !   Lecture du premier etat des reanalyses.         !   Lecture du premier etat des reanalyses.
213         call read_reanalyse(1,ps &         CALL read_reanalyse(1, ps, ucovrea2, vcovrea2, tetarea2, qrea2, &
214              ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)              masserea2, psrea2, 1, nlev)
215         qrea2(:,:)=max(qrea2(:,:),0.1)         qrea2(:, :) = max(qrea2(:, :), 0.1)
216    
217    
        !-----------------------------------------------------------------------  
218         !   Debut de l'integration temporelle:         !   Debut de l'integration temporelle:
219         !   ----------------------------------      END IF ! first
220    
221      endif ! first      ! IMPORTATION DES VENTS, PRESSION ET TEMPERATURE REELS:
     !  
     !-----------------------------------------------------------------------  
     !----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS:  
     !-----------------------------------------------------------------------  
222    
223      ditau=real(itau)      ditau = real(itau)
224      DDAY_step=real(day_step)      dday_step = real(day_step)
225      write(*,*)'ditau,dday_step'      WRITE (*, *) 'ditau, dday_step'
226      write(*,*)ditau,dday_step      WRITE (*, *) ditau, dday_step
227      toto=4*ditau/dday_step      toto = 4*ditau/dday_step
228      reste=toto-aint(toto)      reste = toto - aint(toto)
229      !     write(*,*)'toto,reste',toto,reste  
230        IF (reste==0.) THEN
231      if (reste.eq.0.) then         IF (itau_test==itau) THEN
232         if (itau_test.eq.itau) then            WRITE (*, *) 'deuxieme passage de advreel a itau=', itau
233            write(*,*)'deuxieme passage de advreel a itau=',itau            STOP
234            stop         ELSE
235         else            vcovrea1(:, :) = vcovrea2(:, :)
236            vcovrea1(:,:)=vcovrea2(:,:)            ucovrea1(:, :) = ucovrea2(:, :)
237            ucovrea1(:,:)=ucovrea2(:,:)            tetarea1(:, :) = tetarea2(:, :)
238            tetarea1(:,:)=tetarea2(:,:)            qrea1(:, :) = qrea2(:, :)
239            qrea1(:,:)=qrea2(:,:)  
240              PRINT *, 'LECTURE REANALYSES, pas ', step_rea, 'apres ', &
241            print*,'LECTURE REANALYSES, pas ',step_rea &                 count_no_rea, ' non lectures'
242                 ,'apres ',count_no_rea,' non lectures'            step_rea = step_rea + 1
243            step_rea=step_rea+1            itau_test = itau
244            itau_test=itau            CALL read_reanalyse(step_rea, ps, ucovrea2, vcovrea2, tetarea2, &
245            call read_reanalyse(step_rea,ps &                 qrea2, masserea2, psrea2, 1, nlev)
246                 ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)            qrea2(:, :) = max(qrea2(:, :), 0.1)
247            qrea2(:,:)=max(qrea2(:,:),0.1)            factt = dtvr*iperiod/daysec
248            factt=dtvr*iperiod/daysec            ztau(:) = factt/max(alpha_t(:), 1.E-10)
249            ztau(:)=factt/max(alpha_T(:),1.e-10)            CALL wrgrads(igrads, 1, aire, 'aire      ', 'aire      ')
250            call wrgrads(igrads,1,aire   ,'aire      ','aire      ' )            CALL wrgrads(igrads, 1, dxdys, 'dxdy      ', 'dxdy      ')
251            call wrgrads(igrads,1,dxdys  ,'dxdy      ','dxdy      ' )            CALL wrgrads(igrads, 1, alpha_u, 'au        ', 'au        ')
252            call wrgrads(igrads,1,alpha_u,'au        ','au        ' )            CALL wrgrads(igrads, 1, alpha_t, 'at        ', 'at        ')
253            call wrgrads(igrads,1,alpha_T,'at        ','at        ' )            CALL wrgrads(igrads, 1, ztau, 'taut      ', 'taut      ')
254            call wrgrads(igrads,1,ztau,'taut      ','taut      ' )            CALL wrgrads(igrads, llm, ucov, 'u         ', 'u         ')
255            call wrgrads(igrads,llm,ucov,'u         ','u         ' )            CALL wrgrads(igrads, llm, ucovrea2, 'ua        ', 'ua        ')
256            call wrgrads(igrads,llm,ucovrea2,'ua        ','ua        ' )            CALL wrgrads(igrads, llm, teta, 'T         ', 'T         ')
257            call wrgrads(igrads,llm,teta,'T         ','T         ' )            CALL wrgrads(igrads, llm, tetarea2, 'Ta        ', 'Ta        ')
258            call wrgrads(igrads,llm,tetarea2,'Ta        ','Ta        ' )            CALL wrgrads(igrads, llm, qrea2, 'Qa        ', 'Qa        ')
259            call wrgrads(igrads,llm,qrea2,'Qa        ','Qa        ' )            CALL wrgrads(igrads, llm, q, 'Q         ', 'Q         ')
260            call wrgrads(igrads,llm,q,'Q         ','Q         ' )  
261              CALL wrgrads(igrads, llm, qsat, 'QSAT      ', 'QSAT      ')
262    
263           END IF
264        ELSE
265           count_no_rea = count_no_rea + 1
266        END IF
267    
           call wrgrads(igrads,llm,qsat,'QSAT      ','QSAT      ' )  
   
        endif  
     else  
        count_no_rea=count_no_rea+1  
     endif  
   
     !-----------------------------------------------------------------------  
268      !   Guidage      !   Guidage
269      !    x_gcm = a * x_gcm + (1-a) * x_reanalyses      !    x_gcm = a * x_gcm + (1-a) * x_reanalyses
     !-----------------------------------------------------------------------  
270    
271      if(ini_anal) print*,'ATTENTION !!! ON PART DU GUIDAGE'      IF (ini_anal) PRINT *, 'ATTENTION !!! ON PART DU GUIDAGE'
272    
273      ditau=real(itau)      ditau = real(itau)
274      dday_step=real(day_step)      dday_step = real(day_step)
275    
276    
277      tau=4*ditau/dday_step      tau = 4*ditau/dday_step
278      tau=tau-aint(tau)      tau = tau - aint(tau)
279    
280      !  ucov      !  ucov
281      if (guide_u) then      IF (guide_u) THEN
282         do l=1,llm         DO l = 1, llm
283            do ij=1,ip1jmp1            DO ij = 1, ip1jmp1
284               a=(1.-tau)*ucovrea1(ij,l)+tau*ucovrea2(ij,l)               a = (1.-tau)*ucovrea1(ij, l) + tau*ucovrea2(ij, l)
285               ucov(ij,l)=(1.-alpha_u(ij))*ucov(ij,l)+alpha_u(ij)*a               ucov(ij, l) = (1.-alpha_u(ij))*ucov(ij, l) + alpha_u(ij)*a
286               if (first.and.ini_anal) ucov(ij,l)=a               IF (first .AND. ini_anal) ucov(ij, l) = a
287            enddo            END DO
288         enddo         END DO
289      endif      END IF
290    
291      !  teta      IF (guide_t) THEN
292      if (guide_T) then         DO l = 1, llm
293         do l=1,llm            DO ij = 1, ip1jmp1
294            do ij=1,ip1jmp1               a = (1.-tau)*tetarea1(ij, l) + tau*tetarea2(ij, l)
295               a=(1.-tau)*tetarea1(ij,l)+tau*tetarea2(ij,l)               teta(ij, l) = (1.-alpha_t(ij))*teta(ij, l) + alpha_t(ij)*a
296               teta(ij,l)=(1.-alpha_T(ij))*teta(ij,l)+alpha_T(ij)*a               IF (first .AND. ini_anal) teta(ij, l) = a
297               if (first.and.ini_anal) teta(ij,l)=a            END DO
298            enddo         END DO
299         enddo      END IF
     endif  
300    
301      !  P      !  P
302      if (guide_P) then      IF (guide_p) THEN
303         do ij=1,ip1jmp1         DO ij = 1, ip1jmp1
304            a=(1.-tau)*psrea1(ij)+tau*psrea2(ij)            a = (1.-tau)*psrea1(ij) + tau*psrea2(ij)
305            ps(ij)=(1.-alpha_P(ij))*ps(ij)+alpha_P(ij)*a            ps(ij) = (1.-alpha_p(ij))*ps(ij) + alpha_p(ij)*a
306            if (first.and.ini_anal) ps(ij)=a            IF (first .AND. ini_anal) ps(ij) = a
307         enddo         END DO
308         CALL pression(ip1jmp1,ap,bp,ps,p)         forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
309         CALL massdair(p,masse)         CALL massdair(p, masse)
310      endif      END IF
311    
312    
313      !  q      !  q
314      if (guide_Q) then      IF (guide_q) THEN
315         do l=1,llm         DO l = 1, llm
316            do ij=1,ip1jmp1            DO ij = 1, ip1jmp1
317               a=(1.-tau)*qrea1(ij,l)+tau*qrea2(ij,l)               a = (1.-tau)*qrea1(ij, l) + tau*qrea2(ij, l)
318               !   hum relative en % -> hum specif               !   hum relative en % -> hum specif
319               a=qsat(ij,l)*a*0.01               a = qsat(ij, l)*a*0.01
320               q(ij,l)=(1.-alpha_Q(ij))*q(ij,l)+alpha_Q(ij)*a               q(ij, l) = (1.-alpha_q(ij))*q(ij, l) + alpha_q(ij)*a
321               if (first.and.ini_anal) q(ij,l)=a               IF (first .AND. ini_anal) q(ij, l) = a
322            enddo            END DO
323         enddo         END DO
324      endif      END IF
325    
326      ! vcov      ! vcov
327      if (guide_v) then      IF (guide_v) THEN
328         do l=1,llm         DO l = 1, llm
329            do ij=1,ip1jm            DO ij = 1, ip1jm
330               a=(1.-tau)*vcovrea1(ij,l)+tau*vcovrea2(ij,l)               a = (1.-tau)*vcovrea1(ij, l) + tau*vcovrea2(ij, l)
331               vcov(ij,l)=(1.-alpha_v(ij))*vcov(ij,l)+alpha_v(ij)*a               vcov(ij, l) = (1.-alpha_v(ij))*vcov(ij, l) + alpha_v(ij)*a
332               if (first.and.ini_anal) vcov(ij,l)=a               IF (first .AND. ini_anal) vcov(ij, l) = a
333            enddo            END DO
334            if (first.and.ini_anal) vcov(ij,l)=a            IF (first .AND. ini_anal) vcov(ij, l) = a
335         enddo         END DO
336      endif      END IF
   
     !     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.  
        else  
           gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)  
           print*,'gamma=',gamma  
           if (gamma.lt.1.e-5) then  
              print*,'gamma =',gamma,'<1e-5'  
              stop  
           endif  
           print*,'gamma=',gamma  
           gamma=log(0.5)/log(gamma)  
        endif  
     endif  
   
     alphamin=factt/taumax  
     alphamax=factt/taumin  
   
     do j=1,pjm  
        do i=1,pim  
           if (type.eq.1) then  
              dxdy_=dxdys(i,j)  
              zlat=rlatu(j)*180./pi  
           elseif (type.eq.2) then  
              dxdy_=dxdyu(i,j)  
              zlat=rlatu(j)*180./pi  
           elseif (type.eq.3) then  
              dxdy_=dxdyv(i,j)  
              zlat=rlatv(j)*180./pi  
           endif  
           if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.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.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  
337    
338        first = .FALSE.
339    
340      return    END SUBROUTINE guide
   end subroutine tau2alpha  
341    
342  end module guide_m  END MODULE guide_m

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

  ViewVC Help
Powered by ViewVC 1.1.21