source: LMDZ4/trunk/libf/dyn3d/guide.F @ 802

Last change on this file since 802 was 762, checked in by Laurent Fairhead, 17 years ago

Merge entre la version V3_conv et le HEAD
YM, JG, LF

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