source: LMDZ4/branches/V3_test/libf/dyn3dpar/guide_p.F @ 3925

Last change on this file since 3925 was 709, checked in by Laurent Fairhead, 18 years ago

Nouvelles versions de la dynamique YM
LF

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