source: LMDZ4/branches/IPSL-CM4_IPCC_branch/libf/dyn3d/read_reanalyse.F @ 3046

Last change on this file since 3046 was 524, checked in by lmdzadmin, 21 years ago

Initial revision

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