source: LMDZ.3.3/trunk/libf/dyn3d/read_reanalyse.F @ 604

Last change on this file since 604 was 356, checked in by lmdz, 23 years ago

Lecture du fichier netcdf pour determination du nombre de niveaux verticaux MAF
LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 9.6 KB
Line 
1c
2c $Header$
3c
4      subroutine read_reanalyse(timestep,u,v,t,masse,ps,mode,nlevnc)
5
6c   mode=0 variables naturelles
7c   mode=1 variabels GCM
8
9      IMPLICIT NONE
10#include "netcdf.inc"
11#include "dimensions.h"
12#include "paramet.h"
13#include "comgeom.h"
14#include "comvert.h"
15
16      integer nlevnc
17c      parameter (nlevnc=15)
18C pour annee 2000     
19c      parameter (nlevnc=21)
20      integer timestep,mode,l
21
22      real u(iip1,jjp1,llm),v(iip1,jjm,llm)
23      real t(iip1,jjp1,llm),ps(iip1,jjp1)
24      real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
25
26      integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
27      integer varidpl
28      save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
29      save varidpl
30
31      real*4 unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
32      real*4 tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
33      real*4 pl(nlevnc)
34
35      integer start(4),count(4),status
36
37      real rcode
38      logical first,ncep
39      save first,ncep
40cMod 11 2 99
41       data first,ncep/.true.,.false./
42c      data first,ncep/.true.,.true./
43      print*,'ok0'
44c -----------------------------------------------------------------
45c   Initialisation de la lecture des fichiers
46c -----------------------------------------------------------------
47      if (first) then
48            ncidu=NCOPN('u.nc',NCNOWRIT,rcode)
49            varidu=NCVID(ncidu,'UWND',rcode)
50            print*,'ncidu,varidu',ncidu,varidu
51         if (ncep) then
52            varidpl=NCVID(ncidu,'LEVEL',rcode)
53         else
54            varidpl=NCVID(ncidu,'PRESSURE',rcode)
55         endif
56            print*,'ncidu,varidpl',ncidu,varidpl
57            ncidv=NCOPN('v.nc',NCNOWRIT,rcode)
58            varidv=NCVID(ncidv,'VWND',rcode)
59            print*,'ncidv,varidv',ncidv,varidv
60            ncidt=NCOPN('T.nc',NCNOWRIT,rcode)
61            varidt=NCVID(ncidt,'AIR',rcode)
62            print*,'ncidt,varidt',ncidt,varidt
63c           ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
64c           varidps=NCVID(ncidps,'SP',rcode)
65c           print*,'ncidps,varidps',ncidps,varidps
66      endif
67
68      print*,'ok1'
69c -----------------------------------------------------------------
70c   lecture des champs u, v, T, ps
71c -----------------------------------------------------------------
72
73c  niveaux de pression
74      print*,'WARNING!!! Il n y a pas de test de coherence'
75      print*,'sur le nombre de niveaux verticaux dans le fichier nc'
76#ifdef NC_DOUBLE
77      status=NF_GET_VARA_DOUBLE(ncidu,varidpl,1,nlevnc,pl)
78#else
79      status=NF_GET_VARA_REAL(ncidu,varidpl,1,nlevnc,pl)
80#endif
81      start(1)=1
82      start(2)=1
83      start(3)=1
84      start(4)=timestep
85
86      count(1)=iip1
87      count(2)=jjp1
88      count(3)=nlevnc
89      count(4)=1
90
91c  u
92#ifdef NC_DOUBLE
93      status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unc)
94#else
95      status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)
96#endif
97      print*,'WARNING!!! Correction bidon pour palier a un '
98      print*,'probleme dans la creation des fichiers nc'
99      call correctbid(iim,jjp1*nlevnc,unc)
100c     call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ')
101c  T
102#ifdef NC_DOUBLE
103      status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnc)
104#else
105      status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnc)
106#endif
107      call correctbid(iim,jjp1*nlevnc,tnc)
108c     call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 ')
109
110      count(2)=jjm
111c  v
112#ifdef NC_DOUBLE
113      status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnc)
114#else
115      status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc)
116#endif
117      call correctbid(iim,jjm*nlevnc,vnc)
118c     call dump2d(iip1,jjm,vnc,'VNC COUCHE 1 ')
119
120      start(3)=timestep
121      start(4)=0
122      count(2)=jjp1
123      count(3)=1
124      count(4)=0
125c  ps
126C#ifdef NC_DOUBLE
127c     status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnc)
128C#else
129c     status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)
130C#endif
131c     call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ')
132c     call correctbid(iim,jjp1,psnc)
133
134c   Transformations
135      do l=1,nlevnc
136         pl(l)=pl(l)*100.
137      enddo
138      if(first) then
139         do l=1,nlevnc
140            print*,'PL(',l,')=',pl(l)
141         enddo
142      endif
143
144      call reanalyse2nat(nlevnc,unc,vnc,tnc,psnc,pl,u,v,t
145     s    ,ps,masse,pk)
146c     call dump2d(iip1,jjm,v,'V COUCHE 1 ')
147
148      if(mode.eq.1) call nat2gcm(u,v,t,pk,u,v,t)
149      print*,'TIMESTEP ',timestep
150      if(mode.ne.1) stop'mode pas egal 0'
151c     call dump2d(iip1,jjm,v,'VCOV COUCHE 1 ')
152
153c   Lignes introduites a une epoque pour un probleme oublie...
154c     do l=1,llm
155c        do i=1,iip1
156c           v(i,1,l)=0.
157c           v(i,jjm,l)=0.
158c        enddo
159c     enddo
160      first=.false.
161
162      return
163      end
164
165c -----------------------------------------------------------------
166      subroutine reanalyse2nat(nlevnc,unc,vnc,tnc,psnc,pl,u,v,t
167     s   ,ps,masse,pk)
168
169c -----------------------------------------------------------------
170c   Inversion Nord/sud de la grille + interpollation sur les niveaux
171c   verticaux du modele.
172c -----------------------------------------------------------------
173
174      implicit none
175
176#include "dimensions.h"
177#include "paramet.h"
178#include "comgeom2.h"
179#include "comvert.h"
180#include "comconst.h"
181
182      integer nlevnc
183      real u(iip1,jjp1,llm),v(iip1,jjm,llm)
184      real t(iip1,jjp1,llm),ps(iip1,jjp1)
185
186      real pl(nlevnc)
187      real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
188      real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
189
190      real zu(iip1,jjp1,llm),zv(iip1,jjm,llm)
191      real zt(iip1,jjp1,llm)
192
193      real pext(iip1,jjp1,llm)
194      real pbarx(iip1,jjp1,llm),pbary(iip1,jjm,llm)
195      real plunc(iip1,jjp1,llm),plvnc(iip1,jjm,llm)
196      real plsnc(iip1,jjp1,llm)
197
198      REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
199      real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm),pks(iip1,jjp1)
200      real pkf(iip1,jjp1,llm)
201      real masse(iip1,jjp1,llm),pls(iip1,jjp1,llm)
202      real prefkap,unskap
203
204      integer i,j,l
205
206
207c -----------------------------------------------------------------
208c   calcul de la pression au milieu des couches.
209c -----------------------------------------------------------------
210
211c     do j=1,jjp1
212c        do i=1,iim
213c           ps(i,j)=psnc(i,jjp1+1-j)
214c        enddo
215c        ps(iip1,j)=ps(1,j)
216c     enddo
217
218      CALL pression( ip1jmp1, ap, bp, ps, p )
219      call massdair(p,masse)
220      CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
221c    ....  Calcul de pls , pression au milieu des couches ,en Pascals
222      unskap=1./kappa
223      prefkap =  preff  ** kappa
224c     PRINT *,' Pref kappa unskap  ',preff,kappa,unskap
225      DO l = 1, llm
226       DO j=1,jjp1
227        DO i =1, iip1
228        pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
229        ENDDO
230       ENDDO
231       ENDDO
232
233
234
235c -----------------------------------------------------------------
236c   calcul des pressions pour les grilles u et v
237c -----------------------------------------------------------------
238
239      do l=1,llm
240      do j=1,jjp1
241         do i=1,iip1
242            pext(i,j,l)=pls(i,j,l)*aire(i,j)
243         enddo
244      enddo
245      enddo
246      call massbar(pext, pbarx, pbary )
247      do l=1,llm
248      do j=1,jjp1
249         do i=1,iip1
250            plunc(i,jjp1+1-j,l)=pbarx(i,j,l)/aireu(i,j)
251            plsnc(i,jjp1+1-j,l)=pls(i,j,l)
252         enddo
253      enddo
254      enddo
255      do l=1,llm
256      do j=1,jjm
257         do i=1,iip1
258            plvnc(i,jjm+1-j,l)=pbary(i,j,l)/airev(i,j)
259         enddo
260      enddo
261      enddo
262
263c -----------------------------------------------------------------
264      call pres2lev(unc,zu,nlevnc,llm,pl,plunc,iip1,jjp1)
265      call pres2lev(vnc,zv,nlevnc,llm,pl,plvnc,iip1,jjm )
266      call pres2lev(tnc,zt,nlevnc,llm,pl,plsnc,iip1,jjp1)
267
268c     call dump2d(iip1,jjp1,ps,'PS    ')
269c     call dump2d(iip1,jjp1,psu,'PS    ')
270c     call dump2d(iip1,jjm,psv,'PS    ')
271c  Inversion Nord/Sud
272      do l=1,llm
273         do j=1,jjp1
274            do i=1,iim
275               u(i,j,l)=zu(i,jjp1+1-j,l)
276               t(i,j,l)=zt(i,jjp1+1-j,l)
277            enddo
278            u(iip1,j,l)=u(1,j,l)
279            t(iip1,j,l)=t(1,j,l)
280         enddo
281      enddo
282
283      do l=1,llm
284         do j=1,jjm
285            do i=1,iim
286               v(i,j,l)=zv(i,jjm+1-j,l)
287            enddo
288            v(iip1,j,l)=v(1,j,l)
289         enddo
290      enddo
291
292      return
293      end
294
295      subroutine nat2gcm(u,v,t,pk,ucov,vcov,teta)
296
297      implicit none
298
299#include "dimensions.h"
300#include "paramet.h"
301#include "comgeom2.h"
302#include "comconst.h"
303#include "comvert.h"
304
305      real u(iip1,jjp1,llm),v(iip1,jjm,llm)
306      real t(iip1,jjp1,llm),pk(iip1,jjp1,llm)
307
308      real ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm)
309      real teta(iip1,jjp1,llm)
310
311      integer i,j,l
312
313
314c   calcul de ucov et de la temperature potentielle
315      do l=1,llm
316         do j=1,jjp1
317            do i=1,iim
318               ucov(i,j,l)=u(i,j,l)*cu(i,j)
319               teta(i,j,l)=t(i,j,l)*cpp/pk(i,j,l)
320            enddo
321            ucov(iip1,j,l)=ucov(1,j,l)
322            teta(iip1,j,l)=teta(1,j,l)
323         enddo
324         do i=1,iip1
325            ucov(i,1,l)=0.
326            ucov(i,jjp1,l)=0.
327            teta(i,1,l)=teta(1,1,l)
328            teta(i,jjp1,l)=teta(1,jjp1,l)
329         enddo
330      enddo
331
332c   calcul de ucov
333      do l=1,llm
334         do j=1,jjm
335            do i=1,iim
336               vcov(i,j,l)=v(i,j,l)*cv(i,j)
337            enddo
338            vcov(iip1,j,l)=vcov(1,j,l)
339         enddo
340      enddo
341
342c     call dump2d(iip1,jjp1,teta,'TETA EN BAS   ')
343c     call dump2d(iip1,jjp1,teta(1,1,llm),'TETA EN HAUT   ')
344
345      return
346      end
347      subroutine correctbid(iim,nl,x)
348      integer iim,nl
349      real x(iim+1,nl)
350      integer i,l
351      real zz
352
353      do l=1,nl
354         do i=2,iim-1
355            if(abs(x(i,l)).gt.1.e10) then
356               zz=0.5*(x(i-1,l)+x(i+1,l))
357c              print*,'correction ',i,l,x(i,l),zz
358               x(i,l)=zz
359            endif
360         enddo
361      enddo
362      return
363      end
Note: See TracBrowser for help on using the repository browser.