source: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/read_reanalyse.F @ 1086

Last change on this file since 1086 was 1046, checked in by lmdzadmin, 16 years ago

Modifs guidage pour utiliser des champs de guidage sur niveaux hybrides
Ajout cle logique ok_gradsfile (.false. par defaut) pour activer sorties grads du guidage
FC/IM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.7 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
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,varidap,varidbp
44      save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
45      save ncidpl,varidap,varidbp
46      save varidpl,ncidQ,varidQ
47
48      real*4 unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
49      real*4 tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
50      real*4 Qnc(iip1,jjp1,nlevnc)
51      real*4 plnc(iip1,jjp1,nlevnc)
52      real*4 apnc(nlevnc),bpnc(nlevnc)
53
54      integer start(4),count(4),status
55      integer i,j,l
56
57      real rcode
58      logical first
59      save first
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! endif (first)
131      endif
132      print*,'ok1'
133
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#ifdef NC_DOUBLE
166      status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unc)
167#else
168      status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)
169#endif
170      print*,'WARNING!!! Correction bidon pour palier a un '
171      print*,'probleme dans la creation des fichiers nc'
172      call correctbid(iim,jjp1*nlevnc,unc)
173      call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ')
174      endif
175
176c  Temperature
177c  -----------
178
179      print*,'ncidt=',ncidt,'varidt=',varidt,'start=',start
180      print*,'count=',count
181      if (guide_T) then
182#ifdef NC_DOUBLE
183      status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnc)
184#else
185      status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnc)
186#endif
187      call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 AAA ')
188      call correctbid(iim,jjp1*nlevnc,tnc)
189      call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 BBB ')
190      endif
191
192c  Humidite
193c  --------
194
195      if (guide_Q) then
196#ifdef NC_DOUBLE
197      status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,Qnc)
198#else
199      status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,Qnc)
200#endif
201      call correctbid(iim,jjp1*nlevnc,Qnc)
202      call dump2d(iip1,jjp1,Qnc,'QNC COUCHE 1 ')
203      endif
204
205      count(2)=jjm
206c  Vent meridien
207c  -------------
208
209      if (guide_v) then
210#ifdef NC_DOUBLE
211      status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnc)
212#else
213      status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc)
214#endif
215      call correctbid(iim,jjm*nlevnc,vnc)
216      call dump2d(iip1,jjm,vnc,'VNC COUCHE 1 ')
217      endif
218
219      start(3)=timestep
220      start(4)=0
221      count(2)=jjp1
222      count(3)=1
223      count(4)=0
224
225c  Pression de surface
226c  -------------------
227      psnc(:,:)=0.
228      if ((guide_P).OR.(guide_modele))  then
229#ifdef NC_DOUBLE
230      status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnc)
231#else
232      status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)
233#endif
234      call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ')
235      call correctbid(iim,jjp1,psnc)
236      endif
237
238c Calcul de la pression aux differents niveaux
239c --------------------------------------------
240      if (guide_modele) then
241#ifdef NC_DOUBLE
242      status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
243      status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
244#else
245      status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
246      status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
247#endif
248      else
249#ifdef NC_DOUBLE
250      status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
251#else
252      status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
253#endif
254c conversion en Pascals
255      apnc=apnc*100.
256      bpnc(:)=0.
257      endif
258      do i=1,iip1
259        do j=1,jjp1
260          do l=1,nlevnc
261            plnc(i,j,l)=apnc(l)+bpnc(l)*psnc(i,j)
262          enddo
263        enddo
264      enddo
265      if (first) then
266        print*,'Verification ordre niveaux verticaux'
267        print*,'LMDZ :'
268        do l=1,llm
269           print*,'PL(',l,')=',(ap(l)+ap(l+1))/2.
270     $                    +psi(1,jjp1)*(bp(l)+bp(l+1))/2.
271        enddo
272        print*,'Fichiers guidage'
273        do l=1,nlevnc
274           print*,'PL(',l,')=',plnc(1,1,l)
275        enddo
276        if (guide_u) then
277          do l=1,nlevnc
278            print*,'U(',l,')=',unc(1,1,l)
279          enddo
280        endif
281        if (guide_T) then
282          do l=1,nlevnc
283            print*,'T(',l,')=',tnc(1,1,l)
284          enddo
285        endif
286        print *,'inversion de l''ordre: ok_invertp=',ok_invertp
287      endif
288
289c -----------------------------------------------------------------
290c  Interpollation verticale sur les niveaux modele
291c -----------------------------------------------------------------
292      call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,plnc,u,v,t,Q
293     s    ,ps,masse,pk)
294
295      call dump2d(iip1,jjm,v,'V COUCHE APRES ')
296
297
298c -----------------------------------------------------------------
299c  Passage aux variables du modele (vents covariants, temperature
300c  potentielle et humidite specifique)
301c -----------------------------------------------------------------
302      call nat2gcm(u,v,t,Q,pk,u,v,t,Q)
303      print*,'TIMESTEP ',timestep
304      if(mode.ne.1) stop'mode pas egal 0'
305c     call dump2d(iip1,jjm,v,'VCOV COUCHE 1 ')
306
307c   Lignes introduites a une epoque pour un probleme oublie...
308c     do l=1,llm
309c        do i=1,iip1
310c           v(i,1,l)=0.
311c           v(i,jjm,l)=0.
312c        enddo
313c     enddo
314      first=.false.
315
316      return
317      end
318
319
320c===========================================================================
321      subroutine reanalyse2nat(nlevnc,psi
322     s   ,unc,vnc,tnc,qnc,psnc,plnc,u,v,t,q
323     s   ,ps,masse,pk)
324c===========================================================================
325
326c -----------------------------------------------------------------
327c   Inversion Nord/sud de la grille + interpollation sur les niveaux
328c   verticaux du modele.
329c -----------------------------------------------------------------
330
331      implicit none
332
333#include "dimensions.h"
334#include "paramet.h"
335#include "comgeom2.h"
336#include "comvert.h"
337#include "comconst.h"
338#include "guide.h"
339
340      integer nlevnc
341      real psi(iip1,jjp1)
342      real u(iip1,jjp1,llm),v(iip1,jjm,llm)
343      real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)
344
345      real plnc(iip1,jjp1,nlevnc)
346      real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
347      real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
348      real qnc(iip1,jjp1,nlevnc)
349
350      real zu(iip1,jjp1,llm),zv(iip1,jjm,llm)
351      real zt(iip1,jjp1,llm),zq(iip1,jjp1,llm)
352
353      real pext(iip1,jjp1,llm)
354      real pbarx(iip1,jjp1,llm),pbary(iip1,jjm,llm)
355      real plunc(iip1,jjp1,llm),plvnc(iip1,jjm,llm)
356      real plsnc(iip1,jjp1,llm)
357
358      REAL alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
359      real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm),pks(iip1,jjp1)
360      real pkf(iip1,jjp1,llm)
361      real masse(iip1,jjp1,llm),pls(iip1,jjp1,llm)
362      real prefkap,unskap
363
364
365      integer i,j,l
366
367
368c -----------------------------------------------------------------
369c   calcul de la pression au milieu des couches.
370c -----------------------------------------------------------------
371
372      CALL pression( ip1jmp1, ap, bp, psi, p )
373      call massdair(p,masse)
374      CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
375
376c    ....  Calcul de pls , pression au milieu des couches ,en Pascals
377      unskap=1./kappa
378      prefkap =  preff  ** kappa
379c     PRINT *,' Pref kappa unskap  ',preff,kappa,unskap
380      DO l = 1, llm
381       DO j=1,jjp1
382        DO i =1, iip1
383        pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
384        ENDDO
385       ENDDO
386       ENDDO
387
388
389c -----------------------------------------------------------------
390c   calcul des pressions pour les grilles u et v
391c -----------------------------------------------------------------
392
393      do l=1,llm
394      do j=1,jjp1
395         do i=1,iip1
396            pext(i,j,l)=pls(i,j,l)*aire(i,j)
397         enddo
398      enddo
399      enddo
400      call massbar(pext, pbarx, pbary )
401      do l=1,llm
402      do j=1,jjp1
403         do i=1,iip1
404            plunc(i,jjp1+1-j,l)=pbarx(i,j,l)/aireu(i,j)
405            plsnc(i,jjp1+1-j,l)=pls(i,j,l)
406         enddo
407      enddo
408      enddo
409      do l=1,llm
410      do j=1,jjm
411         do i=1,iip1
412            plvnc(i,jjm+1-j,l)=pbary(i,j,l)/airev(i,j)
413         enddo
414      enddo
415      enddo
416
417c -----------------------------------------------------------------
418
419      if (guide_P) then
420      do j=1,jjp1
421         do i=1,iim
422            ps(i,j)=psnc(i,jjp1+1-j)
423         enddo
424         ps(iip1,j)=ps(1,j)
425      enddo
426      endif
427
428
429c -----------------------------------------------------------------
430      call pres2lev(unc,zu,nlevnc,llm,plnc,plunc,iip1,jjp1,ok_invertp)
431      call pres2lev(vnc,zv,nlevnc,llm,plnc,plvnc,iip1,jjm,ok_invertp )
432      call pres2lev(tnc,zt,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp)
433      call pres2lev(qnc,zq,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp)
434
435c     call dump2d(iip1,jjp1,ps,'PS    ')
436c     call dump2d(iip1,jjp1,psu,'PS    ')
437c     call dump2d(iip1,jjm,psv,'PS    ')
438c  Inversion Nord/Sud
439      do l=1,llm
440         do j=1,jjp1
441            do i=1,iim
442               u(i,j,l)=zu(i,jjp1+1-j,l)
443               t(i,j,l)=zt(i,jjp1+1-j,l)
444               q(i,j,l)=zq(i,jjp1+1-j,l)
445            enddo
446            u(iip1,j,l)=u(1,j,l)
447            t(iip1,j,l)=t(1,j,l)
448            q(iip1,j,l)=q(1,j,l)
449         enddo
450      enddo
451
452      do l=1,llm
453         do j=1,jjm
454            do i=1,iim
455               v(i,j,l)=zv(i,jjm+1-j,l)
456            enddo
457            v(iip1,j,l)=v(1,j,l)
458         enddo
459      enddo
460
461      return
462      end
463
464c===========================================================================
465      subroutine nat2gcm(u,v,t,rh,pk,ucov,vcov,teta,q)
466c===========================================================================
467
468      implicit none
469
470#include "dimensions.h"
471#include "paramet.h"
472#include "comgeom2.h"
473#include "comconst.h"
474#include "comvert.h"
475#include "guide.h"
476
477      real u(iip1,jjp1,llm),v(iip1,jjm,llm)
478      real t(iip1,jjp1,llm),pk(iip1,jjp1,llm),rh(iip1,jjp1,llm)
479      real ps(iip1,jjp1)
480
481      real ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm)
482      real teta(iip1,jjp1,llm),q(iip1,jjp1,llm)
483
484      real pres(iip1,jjp1,llm),qsat(iip1,jjp1,llm)
485
486      real unskap
487
488      integer i,j,l
489
490
491      print*,'Entree dans nat2gcm'
492c    ucov(:,:,:)=0.
493c    do l=1,llm
494c       ucov(:,2:jjm,l)=u(:,2:jjm,l)*cu(:,2:jjm)
495c    enddo
496c    ucov(iip1,:,:)=ucov(1,:,:)
497
498c    teta(:,:,:)=t(:,:,:)*cpp/pk(:,:,:)
499c    teta(iip1,:,:)=teta(1,:,:)
500     
501c   calcul de ucov et de la temperature potentielle
502      do l=1,llm
503         do j=1,jjp1
504            do i=1,iim
505               ucov(i,j,l)=u(i,j,l)*cu(i,j)
506               teta(i,j,l)=t(i,j,l)*cpp/pk(i,j,l)
507            enddo
508            ucov(iip1,j,l)=ucov(1,j,l)
509            teta(iip1,j,l)=teta(1,j,l)
510         enddo
511         do i=1,iip1
512            ucov(i,1,l)=0.
513            ucov(i,jjp1,l)=0.
514            teta(i,1,l)=teta(1,1,l)
515            teta(i,jjp1,l)=teta(1,jjp1,l)
516         enddo
517      enddo
518
519c   calcul de ucov
520      do l=1,llm
521         do j=1,jjm
522            do i=1,iim
523               vcov(i,j,l)=v(i,j,l)*cv(i,j)
524            enddo
525            vcov(iip1,j,l)=vcov(1,j,l)
526         enddo
527      enddo
528
529c     call dump2d(iip1,jjp1,teta,'TETA EN BAS   ')
530c     call dump2d(iip1,jjp1,teta(1,1,llm),'TETA EN HAUT   ')
531
532c  Humidite relative -> specifique
533c  -------------------------------
534      if (guide_hr) then
535c   FINALEMENT ON GUIDE EN HUMIDITE RELATIVE
536      print*,'calcul de unskap'
537      unskap   = 1./ kappa
538      print*,'calcul de pres'
539      pres(:,:,:)=preff*(pk(:,:,:)/cpp)**unskap
540      print*,'calcul de qsat'
541      call q_sat(iip1*jjp1*llm,t,pres,qsat)
542      print*,'calcul de q'
543c   ATTENTION : humidites relatives en %
544      rh(:,:,:)=max(rh(:,:,:)*0.01,1.e-6)
545      q(:,:,:)=qsat(:,:,:)*rh(:,:,:)
546      print*,'calcul de q OK'
547      call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1   ')
548      call dump2d(iip1,jjp1,pres,'PRESSION PREMIERE COUCHE   ')
549      else
550      q(:,:,:)=rh(:,:,:)
551      print*,'calcul de q OK'
552      call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1   ')
553      endif
554
555      return
556      end
557
558
559
560c===========================================================================
561      subroutine correctbid(iim,nl,x)
562c===========================================================================
563      integer iim,nl
564      real x(iim+1,nl)
565      integer i,l
566      real zz
567
568      do l=1,nl
569         do i=2,iim-1
570            if(abs(x(i,l)).gt.1.e10) then
571               zz=0.5*(x(i-1,l)+x(i+1,l))
572c              print*,'correction ',i,l,x(i,l),zz
573               x(i,l)=zz
574            endif
575         enddo
576      enddo
577      return
578      end
Note: See TracBrowser for help on using the repository browser.