source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/read_reanalyse.F @ 406

Last change on this file since 406 was 232, checked in by lmdzadmin, 23 years ago

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