source: LMDZ4/branches/LMDZ4_V2_patch/libf/dyn3d/guide.F @ 5297

Last change on this file since 5297 was 665, checked in by (none), 19 years ago

This commit was manufactured by cvs2svn to create branch
'LMDZ4_V2_patch'.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.6 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            if (first.and.ini_anal) vcov(ij,l)=a
370         enddo
371      endif
372
373c     call dump2d(iip1,jjp1,tetarea1,'TETA REA 1     ')
374c     call dump2d(iip1,jjp1,tetarea2,'TETA REA 2     ')
375c     call dump2d(iip1,jjp1,teta,'TETA           ')
376
377         first=.false.
378
379      return
380      end
381
382c=======================================================================
383      subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)
384c=======================================================================
385
386      implicit none
387
388#include "dimensions.h"
389#include "paramet.h"
390#include "comconst.h"
391#include "comgeom2.h"
392#include "guide.h"
393#include "serre.h"
394
395c   arguments :
396      integer type
397      integer pim,pjm
398      real factt,taumin,taumax,dxdymin,dxdymax
399      real dxdy_,alpha(pim,pjm)
400      real dxdy_min,dxdy_max
401
402c  local :
403      real alphamin,alphamax,gamma,xi
404      save gamma
405      integer i,j,ilon,ilat
406
407      logical first
408      save first
409      data first/.true./
410
411      real cus(iip1,jjp1),cvs(iip1,jjp1)
412      real cuv(iip1,jjm),cvu(iip1,jjp1)
413      real zdx(iip1,jjp1),zdy(iip1,jjp1)
414
415      real zlat
416      real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
417      common/comdxdy/dxdys,dxdyu,dxdyv
418
419      if (first) then
420         do j=2,jjm
421            do i=2,iip1
422               zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
423            enddo
424            zdx(1,j)=zdx(iip1,j)
425         enddo
426         do j=2,jjm
427            do i=1,iip1
428               zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
429            enddo
430         enddo
431         do i=1,iip1
432            zdx(i,1)=zdx(i,2)
433            zdx(i,jjp1)=zdx(i,jjm)
434            zdy(i,1)=zdy(i,2)
435            zdy(i,jjp1)=zdy(i,jjm)
436         enddo
437         do j=1,jjp1
438            do i=1,iip1
439               dxdys(i,j)=sqrt(zdx(i,j)*zdx(i,j)+zdy(i,j)*zdy(i,j))
440            enddo
441         enddo
442         do j=1,jjp1
443            do i=1,iim
444               dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
445            enddo
446            dxdyu(iip1,j)=dxdyu(1,j)
447         enddo
448         do j=1,jjm
449            do i=1,iip1
450               dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
451            enddo
452         enddo
453
454         call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL  ')
455         call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U     ')
456         call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v     ')
457
458c   coordonnees du centre du zoom
459           call coordij(clon,clat,ilon,ilat)
460c   aire de la maille au centre du zoom
461           dxdy_min=dxdys(ilon,ilat)
462c   dxdy maximale de la maille
463           dxdy_max=0.
464           do j=1,jjp1
465              do i=1,iip1
466                 dxdy_max=max(dxdy_max,dxdys(i,j))
467              enddo
468           enddo
469
470         if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
471             print*,'ATTENTION modele peu zoome'
472             print*,'ATTENTION on prend une constante de guidage cste'
473             gamma=0.
474         else
475            gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
476            print*,'gamma=',gamma
477            if (gamma.lt.1.e-5) then
478              print*,'gamma =',gamma,'<1e-5'
479              stop
480            endif
481            print*,'gamma=',gamma
482            gamma=log(0.5)/log(gamma)
483         endif
484      endif
485
486      alphamin=factt/taumax
487      alphamax=factt/taumin
488
489      do j=1,pjm
490         do i=1,pim
491            if (type.eq.1) then
492               dxdy_=dxdys(i,j)
493               zlat=rlatu(j)*180./pi
494            elseif (type.eq.2) then
495               dxdy_=dxdyu(i,j)
496               zlat=rlatu(j)*180./pi
497            elseif (type.eq.3) then
498               dxdy_=dxdyv(i,j)
499               zlat=rlatv(j)*180./pi
500            endif
501         if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
502c  pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
503             alpha(i,j)=alphamin
504         else
505            xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
506            xi=min(xi,1.)
507            if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then
508               alpha(i,j)=xi*alphamin+(1.-xi)*alphamax
509            else
510               alpha(i,j)=0.
511            endif
512         endif
513         enddo
514      enddo
515
516
517      return
518      end
Note: See TracBrowser for help on using the repository browser.