source: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/read_reanalyse.F @ 1121

Last change on this file since 1121 was 1121, checked in by musat, 15 years ago

Bug declarations variables guidage
IM/FH

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.9 KB
Line 
1!
2! $Header$
3!
4c
5c
6      subroutine read_reanalyse(timestep,psi
7     s   ,u,v,t,q,masse,ps,mode,nlevnc)
8
9c   mode=0 variables naturelles
10c   mode=1 variabels GCM
11
12       USE parallel
13c -----------------------------------------------------------------
14c   Declarations
15c -----------------------------------------------------------------
16      IMPLICIT NONE
17
18c common
19c ------
20
21#include "netcdf.inc"
22#include "dimensions.h"
23#include "paramet.h"
24#include "comgeom.h"
25#include "comvert.h"
26#include "guide.h"
27
28c arguments
29c ---------
30      integer nlevnc
31      integer timestep,mode
32
33      real psi(iip1,jjp1)
34      real u(iip1,jjp1,llm),v(iip1,jjm,llm)
35      real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)
36      real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
37
38c local
39c -----
40      integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
41      integer ncidpl
42      integer varidpl,ncidQ,varidQ,varidap,varidbp
43      save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
44      save ncidpl,varidap,varidbp
45      save varidpl,ncidQ,varidQ
46
47      real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
48      real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
49      real Qnc(iip1,jjp1,nlevnc)
50      real plnc(iip1,jjp1,nlevnc)
51      real apnc(nlevnc),bpnc(nlevnc)
52
53      integer start(4),count(4),status
54      integer i,j,l
55
56      real rcode
57      logical first
58      save first
59      INTEGER ierr
60
61      data first/.true./
62
63c -----------------------------------------------------------------
64c   Initialisation de la lecture des fichiers
65c -----------------------------------------------------------------
66      if (first) then
67           ncidpl=-99
68           print*,'Intitialisation de read reanalsye'
69
70c Niveaux de pression si non constants
71            if (guide_modele) then
72            print *,'Vous êtes entrain de lire des données sur
73     .               niveaux modèle'
74            ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcode)
75            varidap=NCVID(ncidpl,'AP',rcode)
76            varidbp=NCVID(ncidpl,'BP',rcode)
77            print*,'ncidpl,varidap',ncidpl,varidap
78            endif
79
80c Vent zonal
81            if (guide_u) then
82               ncidu=NCOPN('u.nc',NCNOWRIT,rcode)
83               varidu=NCVID(ncidu,'UWND',rcode)
84               print*,'ncidu,varidu',ncidu,varidu
85               if (ncidpl.eq.-99) ncidpl=ncidu
86            endif
87
88c Vent meridien
89            if (guide_v) then
90               ncidv=NCOPN('v.nc',NCNOWRIT,rcode)
91               varidv=NCVID(ncidv,'VWND',rcode)
92               print*,'ncidv,varidv',ncidv,varidv
93               if (ncidpl.eq.-99) ncidpl=ncidv
94            endif
95
96c Temperature
97            if (guide_T) then
98               ncidt=NCOPN('T.nc',NCNOWRIT,rcode)
99               varidt=NCVID(ncidt,'AIR',rcode)
100               print*,'ncidt,varidt',ncidt,varidt
101               if (ncidpl.eq.-99) ncidpl=ncidt
102            endif
103
104c Humidite
105            if (guide_Q) then
106               ncidQ=NCOPN('hur.nc',NCNOWRIT,rcode)
107               varidQ=NCVID(ncidQ,'RH',rcode)
108               print*,'ncidQ,varidQ',ncidQ,varidQ
109               if (ncidpl.eq.-99) ncidpl=ncidQ
110            endif
111
112c Pression de surface
113            if ((guide_P).OR.(guide_modele)) then
114               ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
115               varidps=NCVID(ncidps,'SP',rcode)
116               print*,'ncidps,varidps',ncidps,varidps
117            endif
118
119c Coordonnee verticale
120            if (.not.guide_modele) then
121              if (ncep) then
122               print*,'Vous etes entrain de lire des donnees NCEP'
123               varidpl=NCVID(ncidpl,'LEVEL',rcode)
124            else
125               print*,'Vous etes entrain de lire des donnees ECMWF'
126               varidpl=NCVID(ncidpl,'PRESSURE',rcode)
127              endif
128              print*,'ncidpl,varidpl',ncidpl,varidpl
129            endif
130            print*,'ncidu,varidpl',ncidu,varidpl
131
132      endif
133      print*,'ok1'
134
135c -----------------------------------------------------------------
136c   lecture des champs u, v, T, ps
137c -----------------------------------------------------------------
138
139c  dimensions pour les champs scalaires et le vent zonal
140c  -----------------------------------------------------
141
142      start(1)=1
143      start(2)=1
144      start(3)=1
145      start(4)=timestep
146
147      count(1)=iip1
148      count(2)=jjp1
149      count(3)=nlevnc
150      count(4)=1
151
152c mise a zero des tableaux
153c ------------------------
154       unc(:,:,:)=0.
155       vnc(:,:,:)=0.
156       tnc(:,:,:)=0.
157       Qnc(:,:,:)=0.
158       plnc(:,:,:)=0.
159
160c  Vent zonal
161c  ----------
162
163      if (guide_u) then
164         print*,'avant la lecture de UNCEP nd de niv:',nlevnc
165
166#ifdef NC_DOUBLE
167         status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unc)
168#else
169         status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)
170#endif
171      print*,'WARNING!!! Correction bidon pour palier a un '
172      print*,'probleme dans la creation des fichiers nc'
173      call correctbid(iim,jjp1*nlevnc,unc)
174      call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ')
175      endif
176
177c  Temperature
178c  -----------
179
180      print*,'ncidt=',ncidt,'varidt=',varidt,'start=',start
181      print*,'count=',count
182      if (guide_T) then
183#ifdef NC_DOUBLE
184         status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnc)
185#else
186         status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnc)
187#endif
188         call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 AAA ')
189         call correctbid(iim,jjp1*nlevnc,tnc)
190         call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 BBB ')
191      endif
192
193c  Humidite
194c  --------
195
196      if (guide_Q) then
197#ifdef NC_DOUBLE
198         status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,Qnc)
199#else
200         status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,Qnc)
201#endif
202         call correctbid(iim,jjp1*nlevnc,Qnc)
203         call dump2d(iip1,jjp1,Qnc,'QNC COUCHE 1 ')
204      endif
205
206      count(2)=jjm
207c  Vent meridien
208c  -------------
209
210      if (guide_v) then
211#ifdef NC_DOUBLE
212         status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnc)
213#else
214         status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc)
215#endif
216         call correctbid(iim,jjm*nlevnc,vnc)
217         call dump2d(iip1,jjm,vnc,'VNC COUCHE 1 ')
218      endif
219
220      start(3)=timestep
221      start(4)=0
222      count(2)=jjp1
223      count(3)=1
224      count(4)=0
225
226c  Pression de surface
227c  -------------------
228      psnc(:,:)=0.
229      if ((guide_P).OR.(guide_modele))  then
230#ifdef NC_DOUBLE
231         status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnc)
232#else
233         status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)
234#endif
235         call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ')
236         call correctbid(iim,jjp1,psnc)
237      endif
238
239c Calcul de la pression aux differents niveaux
240c --------------------------------------------
241      if (guide_modele) then
242#ifdef NC_DOUBLE
243      status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
244      status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
245#else
246      status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
247      status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
248#endif
249      else
250#ifdef NC_DOUBLE
251      status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
252#else
253      status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
254#endif
255c conversion en Pascals
256      apnc=apnc*100.
257      bpnc(:)=0.
258      endif
259      do i=1,iip1
260        do j=1,jjp1
261          do l=1,nlevnc
262            plnc(i,j,l)=apnc(l)+bpnc(l)*psnc(i,j)
263          enddo
264        enddo
265      enddo
266      if (first) then
267        print*,'Verification ordre niveaux verticaux'
268        print*,'LMDZ :'
269        do l=1,llm
270           print*,'PL(',l,')=',(ap(l)+ap(l+1))/2.
271     $                    +psi(1,jjp1)*(bp(l)+bp(l+1))/2.
272        enddo
273        print*,'Fichiers guidage'
274        do l=1,nlevnc
275           print*,'PL(',l,')=',plnc(1,1,l)
276        enddo
277        if (guide_u) then
278          do l=1,nlevnc
279            print*,'U(',l,')=',unc(1,1,l)
280          enddo
281        endif
282        if (guide_T) then
283          do l=1,nlevnc
284            print*,'T(',l,')=',tnc(1,1,l)
285          enddo
286        endif
287        print *,'inversion de l''ordre: ok_invertp=',ok_invertp
288      endif
289
290c -----------------------------------------------------------------
291c  Interpollation verticale sur les niveaux modele
292c -----------------------------------------------------------------
293      call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,plnc,u,v,t,Q
294     s    ,ps,masse,pk)
295
296      call dump2d(iip1,jjm,v,'V COUCHE APRES ')
297
298
299c -----------------------------------------------------------------
300c  Passage aux variables du modele (vents covariants, temperature
301c  potentielle et humidite specifique)
302c -----------------------------------------------------------------
303      call nat2gcm(u,v,t,Q,pk,u,v,t,Q)
304      print*,'TIMESTEP ',timestep
305      if(mode.ne.1) stop'mode pas egal 0'
306c     call dump2d(iip1,jjm,v,'VCOV COUCHE 1 ')
307
308c   Lignes introduites a une epoque pour un probleme oublie...
309c     do l=1,llm
310c        do i=1,iip1
311c           v(i,1,l)=0.
312c           v(i,jjm,l)=0.
313c        enddo
314c     enddo
315      first=.false.
316
317      return
318      end
319
320
321c===========================================================================
322      subroutine reanalyse2nat(nlevnc,psi
323     s   ,unc,vnc,tnc,qnc,psnc,plnc,u,v,t,q
324     s   ,ps,masse,pk)
325c===========================================================================
326
327c -----------------------------------------------------------------
328c   Inversion Nord/sud de la grille + interpollation sur les niveaux
329c   verticaux du modele.
330c -----------------------------------------------------------------
331
332      implicit none
333
334#include "dimensions.h"
335#include "paramet.h"
336#include "comgeom2.h"
337#include "comvert.h"
338#include "comconst.h"
339#include "guide.h"
340
341      integer nlevnc
342      real psi(iip1,jjp1)
343      real u(iip1,jjp1,llm),v(iip1,jjm,llm)
344      real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)
345
346      real plnc(iip1,jjp1,nlevnc)
347      real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
348      real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
349      real qnc(iip1,jjp1,nlevnc)
350
351      real zu(iip1,jjp1,llm),zv(iip1,jjm,llm)
352      real zt(iip1,jjp1,llm),zq(iip1,jjp1,llm)
353
354      real pext(iip1,jjp1,llm)
355      real pbarx(iip1,jjp1,llm),pbary(iip1,jjm,llm)
356      real plunc(iip1,jjp1,llm),plvnc(iip1,jjm,llm)
357      real plsnc(iip1,jjp1,llm)
358
359      REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
360      real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm),pks(iip1,jjp1)
361      real pkf(iip1,jjp1,llm)
362      real masse(iip1,jjp1,llm),pls(iip1,jjp1,llm)
363      real prefkap,unskap
364
365
366      integer i,j,l
367
368
369c -----------------------------------------------------------------
370c   calcul de la pression au milieu des couches.
371c -----------------------------------------------------------------
372
373      CALL pression( ip1jmp1, ap, bp, psi, p )
374      call massdair(p,masse)
375      CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
376
377c    ....  Calcul de pls , pression au milieu des couches ,en Pascals
378      unskap=1./kappa
379      prefkap =  preff  ** kappa
380c     PRINT *,' Pref kappa unskap  ',preff,kappa,unskap
381      DO l = 1, llm
382       DO j=1,jjp1
383        DO i =1, iip1
384        pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
385        ENDDO
386       ENDDO
387       ENDDO
388
389
390c -----------------------------------------------------------------
391c   calcul des pressions pour les grilles u et v
392c -----------------------------------------------------------------
393
394      do l=1,llm
395      do j=1,jjp1
396         do i=1,iip1
397            pext(i,j,l)=pls(i,j,l)*aire(i,j)
398         enddo
399      enddo
400      enddo
401      call massbar(pext, pbarx, pbary )
402      do l=1,llm
403      do j=1,jjp1
404         do i=1,iip1
405            plunc(i,jjp1+1-j,l)=pbarx(i,j,l)/aireu(i,j)
406            plsnc(i,jjp1+1-j,l)=pls(i,j,l)
407         enddo
408      enddo
409      enddo
410      do l=1,llm
411      do j=1,jjm
412         do i=1,iip1
413            plvnc(i,jjm+1-j,l)=pbary(i,j,l)/airev(i,j)
414         enddo
415      enddo
416      enddo
417
418c -----------------------------------------------------------------
419
420      if (guide_P) then
421      do j=1,jjp1
422         do i=1,iim
423            ps(i,j)=psnc(i,jjp1+1-j)
424         enddo
425         ps(iip1,j)=ps(1,j)
426      enddo
427      endif
428
429
430c -----------------------------------------------------------------
431      call pres2lev(unc,zu,nlevnc,llm,plnc,plunc,iip1,jjp1,ok_invertp)
432      call pres2lev(vnc,zv,nlevnc,llm,plnc,plvnc,iip1,jjm,ok_invertp )
433      call pres2lev(tnc,zt,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp)
434      call pres2lev(qnc,zq,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp)
435
436c     call dump2d(iip1,jjp1,ps,'PS    ')
437c     call dump2d(iip1,jjp1,psu,'PS    ')
438c     call dump2d(iip1,jjm,psv,'PS    ')
439c  Inversion Nord/Sud
440      do l=1,llm
441         do j=1,jjp1
442            do i=1,iim
443               u(i,j,l)=zu(i,jjp1+1-j,l)
444               t(i,j,l)=zt(i,jjp1+1-j,l)
445               q(i,j,l)=zq(i,jjp1+1-j,l)
446            enddo
447            u(iip1,j,l)=u(1,j,l)
448            t(iip1,j,l)=t(1,j,l)
449            q(iip1,j,l)=q(1,j,l)
450         enddo
451      enddo
452
453      do l=1,llm
454         do j=1,jjm
455            do i=1,iim
456               v(i,j,l)=zv(i,jjm+1-j,l)
457            enddo
458            v(iip1,j,l)=v(1,j,l)
459         enddo
460      enddo
461
462      return
463      end
464
465c===========================================================================
466      subroutine nat2gcm(u,v,t,rh,pk,ucov,vcov,teta,q)
467c===========================================================================
468
469      implicit none
470
471#include "dimensions.h"
472#include "paramet.h"
473#include "comgeom2.h"
474#include "comconst.h"
475#include "comvert.h"
476#include "guide.h"
477
478      real u(iip1,jjp1,llm),v(iip1,jjm,llm)
479      real t(iip1,jjp1,llm),pk(iip1,jjp1,llm),rh(iip1,jjp1,llm)
480      real ps(iip1,jjp1)
481
482      real ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm)
483      real teta(iip1,jjp1,llm),q(iip1,jjp1,llm)
484
485      real pres(iip1,jjp1,llm),qsat(iip1,jjp1,llm)
486
487      real unskap
488
489      integer i,j,l
490
491
492      print*,'Entree dans nat2gcm'
493c    ucov(:,:,:)=0.
494c    do l=1,llm
495c       ucov(:,2:jjm,l)=u(:,2:jjm,l)*cu(:,2:jjm)
496c    enddo
497c    ucov(iip1,:,:)=ucov(1,:,:)
498
499c    teta(:,:,:)=t(:,:,:)*cpp/pk(:,:,:)
500c    teta(iip1,:,:)=teta(1,:,:)
501     
502c   calcul de ucov et de la temperature potentielle
503      do l=1,llm
504         do j=1,jjp1
505            do i=1,iim
506               ucov(i,j,l)=u(i,j,l)*cu(i,j)
507               teta(i,j,l)=t(i,j,l)*cpp/pk(i,j,l)
508            enddo
509            ucov(iip1,j,l)=ucov(1,j,l)
510            teta(iip1,j,l)=teta(1,j,l)
511         enddo
512         do i=1,iip1
513            ucov(i,1,l)=0.
514            ucov(i,jjp1,l)=0.
515            teta(i,1,l)=teta(1,1,l)
516            teta(i,jjp1,l)=teta(1,jjp1,l)
517         enddo
518      enddo
519
520c   calcul de ucov
521      do l=1,llm
522         do j=1,jjm
523            do i=1,iim
524               vcov(i,j,l)=v(i,j,l)*cv(i,j)
525            enddo
526            vcov(iip1,j,l)=vcov(1,j,l)
527         enddo
528      enddo
529
530c     call dump2d(iip1,jjp1,teta,'TETA EN BAS   ')
531c     call dump2d(iip1,jjp1,teta(1,1,llm),'TETA EN HAUT   ')
532
533c  Humidite relative -> specifique
534c  -------------------------------
535      if (guide_hr) then
536c   FINALEMENT ON GUIDE EN HUMIDITE RELATIVE
537      print*,'calcul de unskap'
538      unskap   = 1./ kappa
539      print*,'calcul de pres'
540      pres(:,:,:)=preff*(pk(:,:,:)/cpp)**unskap
541      print*,'calcul de qsat'
542      call q_sat(iip1*jjp1*llm,t,pres,qsat)
543      print*,'calcul de q'
544c   ATTENTION : humidites relatives en %
545      rh(:,:,:)=max(rh(:,:,:)*0.01,1.e-6)
546      q(:,:,:)=qsat(:,:,:)*rh(:,:,:)
547      print*,'calcul de q OK'
548      call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1   ')
549      call dump2d(iip1,jjp1,pres,'PRESSION PREMIERE COUCHE   ')
550      else
551      q(:,:,:)=rh(:,:,:)
552      print*,'calcul de q OK'
553      call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1   ')
554      endif
555
556      return
557      end
558
559
560
561c===========================================================================
562      subroutine correctbid(iim,nl,x)
563c===========================================================================
564      integer iim,nl
565      real x(iim+1,nl)
566      integer i,l
567      real zz
568
569      do l=1,nl
570         do i=2,iim-1
571            if(abs(x(i,l)).gt.1.e10) then
572               zz=0.5*(x(i-1,l)+x(i+1,l))
573c              print*,'correction ',i,l,x(i,l),zz
574               x(i,l)=zz
575            endif
576         enddo
577      enddo
578      return
579      end
Note: See TracBrowser for help on using the repository browser.