source: LMDZ4/trunk/libf/dyn3dpar/guide_p.F @ 1080

Last change on this file since 1080 was 1049, checked in by lmdzadmin, 16 years ago

Petit oubli endif
AC/IM

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