source: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/guide.F @ 1167

Last change on this file since 1167 was 1135, checked in by lguez, 16 years ago

Translated calls using NetCDF 2.4 interface to calls using NetCDF 3.6
Fortran 90 interface.

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