/[lmdze]/trunk/Sources/dyn3d/Guide/Read_reanalyse/read_reanalyse.f
ViewVC logotype

Diff of /trunk/Sources/dyn3d/Guide/Read_reanalyse/read_reanalyse.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/dyn3d/read_reanalyse.f revision 3 by guez, Wed Feb 27 13:16:39 2008 UTC trunk/Sources/dyn3d/Guide/Read_reanalyse/read_reanalyse.f revision 134 by guez, Wed Apr 29 15:47:56 2015 UTC
# Line 1  Line 1 
1  !  module read_reanalyse_m
2  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/read_reanalyse.F,v 1.3 2005/04/15 12:31:21 lmdzadmin Exp $  
3  !    IMPLICIT NONE
4  c  
5  c  contains
6        subroutine read_reanalyse(timestep,psi  
7       s   ,u,v,t,q,masse,ps,mode,nlevnc)    subroutine read_reanalyse(timestep, psi, u, v, t, q, masse, nlevnc)
8    
9  c   mode=0 variables naturelles      ! From LMDZ4/libf/dyn3d/read_reanalyse.F, version 1.3, 2005/04/15 12:31:21
10  c   mode=1 variabels GCM  
11        USE conf_guide_m, ONLY: guide_q, guide_t, guide_u, guide_v, ncep
12  c -----------------------------------------------------------------      USE dimens_m, ONLY: iim, jjm, llm
13  c   Declarations      USE netcdf, ONLY: nf90_get_var, nf90_inq_varid, nf90_nowrite, nf90_open
14  c -----------------------------------------------------------------      USE paramet_m, ONLY: iip1, jjp1
15        use dimens_m      use reanalyse2nat_m, only: reanalyse2nat
16        use paramet_m  
17        use comvert      integer timestep
18        use comgeom      real, intent(in):: psi(iip1, jjp1)
19        use guide_m      real u(iip1, jjp1, llm), v(iip1, jjm, llm)
20        IMPLICIT NONE      real t(iip1, jjp1, llm), q(iip1, jjp1, llm)
21        real masse(iip1, jjp1, llm)
22  c common      integer nlevnc
23  c ------  
24        ! Local:
25        include "netcdf.inc"  
26        integer l
27        real pk(iip1, jjp1, llm)
28  c arguments      integer, save:: ncidu, varidu, ncidv, varidv, ncidt, varidt
29  c ---------      integer, save:: ncidpl
30        integer nlevnc      integer, save:: varidpl, ncidQ, varidQ
31        integer timestep,mode,l      real unc(iip1, jjp1, nlevnc), vnc(iip1, jjm, nlevnc)
32        real tnc(iip1, jjp1, nlevnc)
33        real psi(iip1,jjp1)      real Qnc(iip1, jjp1, nlevnc)
34        real u(iip1,jjp1,llm),v(iip1,jjm,llm)      real pl(nlevnc)
35        real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)      integer start(4), count(4), status
36        real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)      real rcode
37        logical:: first = .true.
38    
39  c local      ! -----------------------------------------------------------------
40  c -----  
41        integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps      !   Initialisation de la lecture des fichiers
42        integer ncidpl  
43        integer varidpl,ncidQ,varidQ      if (first) then
44        save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps         ncidpl=-99
45        save ncidpl         print *, 'Intitialisation de read reanalsye'
46        save varidpl,ncidQ,varidQ  
47           ! Vent zonal
48        real*4 unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)         if (guide_u) then
49        real*4 tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)            rcode=nf90_open('u.nc', nf90_nowrite, ncidu)
50        real*4 Qnc(iip1,jjp1,nlevnc)            rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
51        real*4 pl(nlevnc)            if (ncidpl.eq.-99) ncidpl=ncidu
52           endif
53        integer start(4),count(4),status  
54           ! Vent meridien
55        real rcode         if (guide_v) then
56        logical first            rcode=nf90_open('v.nc', nf90_nowrite, ncidv)
57        save first            rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
58              if (ncidpl.eq.-99) ncidpl=ncidv
59        data first/.true./         endif
60    
61           ! Temperature
62           if (guide_T) then
63  c -----------------------------------------------------------------            rcode=nf90_open('T.nc', nf90_nowrite, ncidt)
64  c   Initialisation de la lecture des fichiers            rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
65  c -----------------------------------------------------------------            if (ncidpl.eq.-99) ncidpl=ncidt
66        if (first) then         endif
67             ncidpl=-99  
68             print*,'Intitialisation de read reanalsye'         ! Humidite
69           if (guide_Q) then
70  c Vent zonal            rcode=nf90_open('hur.nc', nf90_nowrite, ncidQ)
71              if (guide_u) then            rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
72              ncidu=NCOPN('u.nc',NCNOWRIT,rcode)            if (ncidpl.eq.-99) ncidpl=ncidQ
73              varidu=NCVID(ncidu,'UWND',rcode)         endif
74              print*,'ncidu,varidu',ncidu,varidu  
75              if (ncidpl.eq.-99) ncidpl=ncidu         ! Coordonnee verticale
76              endif         if (ncep) then
77              print *, 'Vous etes entrain de lire des donnees NCEP'
78  c Vent meridien            rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
79              if (guide_v) then         else
80              ncidv=NCOPN('v.nc',NCNOWRIT,rcode)            print *, 'Vous etes entrain de lire des donnees ECMWF'
81              varidv=NCVID(ncidv,'VWND',rcode)            rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
82              print*,'ncidv,varidv',ncidv,varidv         endif
83              if (ncidpl.eq.-99) ncidpl=ncidv      endif
84              endif  
85        ! Niveaux de pression
86  c Temperature  
87              if (guide_T) then      ! Warning: il n y a pas de test de coherence sur le nombre de
88              ncidt=NCOPN('T.nc',NCNOWRIT,rcode)      ! niveaux verticaux dans le fichier nc'
89              varidt=NCVID(ncidt,'AIR',rcode)      status=NF90_GET_VAR(ncidpl, varidpl, pl)
90              print*,'ncidt,varidt',ncidt,varidt      !  passage en pascal
91              if (ncidpl.eq.-99) ncidpl=ncidt      pl(:)=100.*pl(:)
92              endif      if (first) then
93           do l=1, nlevnc
94  c Humidite            print *, 'PL(', l, ')=', pl(l)
             if (guide_Q) then  
             ncidQ=NCOPN('hur.nc',NCNOWRIT,rcode)  
             varidQ=NCVID(ncidQ,'RH',rcode)  
             print*,'ncidQ,varidQ',ncidQ,varidQ  
             if (ncidpl.eq.-99) ncidpl=ncidQ  
             endif  
   
 c Pression de surface  
             if (guide_P) then  
             ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)  
             varidps=NCVID(ncidps,'SP',rcode)  
             print*,'ncidps,varidps',ncidps,varidps  
             endif  
   
 c Coordonnee verticale  
             if (ncep) then  
                print*,'Vous etes entrain de lire des donnees NCEP'  
                varidpl=NCVID(ncidpl,'LEVEL',rcode)  
             else  
                print*,'Vous etes entrain de lire des donnees ECMWF'  
                varidpl=NCVID(ncidpl,'PRESSURE',rcode)  
             endif  
             print*,'ncidu,varidpl',ncidu,varidpl  
       endif  
       print*,'ok1'  
   
 c Niveaux de pression  
       print*,'WARNING!!! Il n y a pas de test de coherence'  
       print*,'sur le nombre de niveaux verticaux dans le fichier nc'  
       status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,pl)  
 c  passage en pascal  
       pl(:)=100.*pl(:)  
       if (first) then  
        do l=1,nlevnc  
           print*,'PL(',l,')=',pl(l)  
95         enddo         enddo
96        endif      endif
97    
98        !   lecture des champs u, v, T
99    
100        !  dimensions pour les champs scalaires et le vent zonal
101    
102        start(1)=1
103        start(2)=1
104        start(3)=1
105        start(4)=timestep
106    
107        count(1)=iip1
108        count(2)=jjp1
109        count(3)=nlevnc
110        count(4)=1
111    
112        ! mise a zero des tableaux
113    
114        unc(:, :, :)=0.
115        vnc(:, :, :)=0.
116        tnc(:, :, :)=0.
117        Qnc(:, :, :)=0.
118    
119        !  Vent zonal
120    
121        if (guide_u) then
122           status=NF90_GET_VAR(ncidu, varidu, unc, start, count)
123           ! Warning Correction bidon pour palier a un probleme dans la
124           ! creation des fichiers nc
125           call correctbid(iim, jjp1*nlevnc, unc)
126        endif
127    
128        !  Temperature
129    
130        if (guide_T) then
131           status=NF90_GET_VAR(ncidt, varidt, tnc, start, count)
132           call correctbid(iim, jjp1*nlevnc, tnc)
133        endif
134    
135        !  Humidite
136    
137        if (guide_Q) then
138           status=NF90_GET_VAR(ncidQ, varidQ, Qnc, start, count)
139           call correctbid(iim, jjp1*nlevnc, Qnc)
140        endif
141    
142        count(2)=jjm
143        !  Vent meridien
144    
145        if (guide_v) then
146           status=NF90_GET_VAR(ncidv, varidv, vnc, start, count)
147           call correctbid(iim, jjm*nlevnc, vnc)
148        endif
149    
150        start(3)=timestep
151        start(4)=0
152        count(2)=jjp1
153        count(3)=1
154        count(4)=0
155    
156        !  Interpolation verticale sur les niveaux modele
157    
158        call reanalyse2nat(nlevnc, psi, unc, vnc, tnc, Qnc, pl, u, v, t, Q, &
159             masse, pk)
160    
161        !  Passage aux variables du modele (vents covariants, temperature
162        !  potentielle et humidite specifique)
163    
164        call nat2gcm(u, v, t, Q, pk, u, v, t, Q)
165        first=.false.
166    
167      end subroutine read_reanalyse
168    
169  c -----------------------------------------------------------------  end module read_reanalyse_m
 c   lecture des champs u, v, T, ps  
 c -----------------------------------------------------------------  
   
 c  dimensions pour les champs scalaires et le vent zonal  
 c  -----------------------------------------------------  
   
       start(1)=1  
       start(2)=1  
       start(3)=1  
       start(4)=timestep  
   
       count(1)=iip1  
       count(2)=jjp1  
       count(3)=nlevnc  
       count(4)=1  
   
 c mise a zero des tableaux  
 c ------------------------  
        unc(:,:,:)=0.  
        vnc(:,:,:)=0.  
        tnc(:,:,:)=0.  
        Qnc(:,:,:)=0.  
   
 c  Vent zonal  
 c  ----------  
   
       if (guide_u) then  
       print*,'avant la lecture de UNCEP nd de niv:',nlevnc  
       status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)  
 c     call dump2d(iip1,jjp1,unc,'VENT NCEP   ')  
 c     call dump2d(iip1,40,unc(1,1,nlevnc),'VENT NCEP   ')  
       print*,'WARNING!!! Correction bidon pour palier a un '  
       print*,'probleme dans la creation des fichiers nc'  
       call correctbid(iim,jjp1*nlevnc,unc)  
       call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ')  
       endif  
   
 c  Temperature  
 c  -----------  
   
       print*,'ncidt=',ncidt,'varidt=',varidt,'start=',start  
       print*,'count=',count  
       if (guide_T) then  
       status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnc)  
       call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 AAA ')  
       call correctbid(iim,jjp1*nlevnc,tnc)  
       call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 BBB ')  
       endif  
   
 c  Humidite  
 c  --------  
   
       if (guide_Q) then  
       status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,Qnc)  
       call correctbid(iim,jjp1*nlevnc,Qnc)  
       call dump2d(iip1,jjp1,Qnc,'QNC COUCHE 1 ')  
       endif  
   
       count(2)=jjm  
 c  Vent meridien  
 c  -------------  
   
       if (guide_v) then  
       status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc)  
       call correctbid(iim,jjm*nlevnc,vnc)  
       call dump2d(iip1,jjm,vnc,'VNC COUCHE 1 ')  
       endif  
   
       start(3)=timestep  
       start(4)=0  
       count(2)=jjp1  
       count(3)=1  
       count(4)=0  
   
 c  Pression de surface  
 c  -------------------  
   
       if (guide_P) then  
       status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)  
       call dump2d(iip1,jjp1,psnc,'PSNC COUCHE 1 ')  
       call correctbid(iim,jjp1,psnc)  
       endif  
   
   
   
 c -----------------------------------------------------------------  
 c  Interpollation verticale sur les niveaux modele  
 c -----------------------------------------------------------------  
       call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,pl,u,v,t,Q  
      s    ,ps,masse,pk)  
   
       call dump2d(iip1,jjm,v,'V COUCHE APRES ')  
   
   
 c -----------------------------------------------------------------  
 c  Passage aux variables du modele (vents covariants, temperature  
 c  potentielle et humidite specifique)  
 c -----------------------------------------------------------------  
       call nat2gcm(u,v,t,Q,pk,u,v,t,Q)  
       print*,'TIMESTEP ',timestep  
       if(mode.ne.1) stop'mode pas egal 0'  
 c     call dump2d(iip1,jjm,v,'VCOV COUCHE 1 ')  
   
 c   Lignes introduites a une epoque pour un probleme oublie...  
 c     do l=1,llm  
 c        do i=1,iip1  
 c           v(i,1,l)=0.  
 c           v(i,jjm,l)=0.  
 c        enddo  
 c     enddo  
       first=.false.  
   
       return  
       end  
   
   
 c===========================================================================  
       subroutine reanalyse2nat(nlevnc,psi  
      s   ,unc,vnc,tnc,qnc,psnc,pl,u,v,t,q  
      s   ,ps,masse,pk)  
 c===========================================================================  
   
 c -----------------------------------------------------------------  
 c   Inversion Nord/sud de la grille + interpollation sur les niveaux  
 c   verticaux du modele.  
 c -----------------------------------------------------------------  
   
       use dimens_m  
       use paramet_m  
       use comconst  
       use comvert  
       use comgeom  
       use exner_hyb_m, only: exner_hyb  
       use guide_m  
       use pression_m, only: pression  
   
       implicit none  
   
   
       integer nlevnc  
       real psi(iip1,jjp1)  
       real u(iip1,jjp1,llm),v(iip1,jjm,llm)  
       real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)  
   
       real pl(nlevnc)  
       real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)  
       real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)  
       real qnc(iip1,jjp1,nlevnc)  
   
       real zu(iip1,jjp1,llm),zv(iip1,jjm,llm)  
       real zt(iip1,jjp1,llm),zq(iip1,jjp1,llm)  
   
       real pext(iip1,jjp1,llm)  
       real pbarx(iip1,jjp1,llm),pbary(iip1,jjm,llm)  
       real plunc(iip1,jjp1,llm),plvnc(iip1,jjm,llm)  
       real plsnc(iip1,jjp1,llm)  
   
       real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm),pks(iip1,jjp1)  
       real pkf(iip1,jjp1,llm)  
       real masse(iip1,jjp1,llm),pls(iip1,jjp1,llm)  
       real prefkap,unskap  
   
   
       integer i,j,l  
   
   
 c -----------------------------------------------------------------  
 c   calcul de la pression au milieu des couches.  
 c -----------------------------------------------------------------  
   
       CALL pression( ip1jmp1, ap, bp, psi, p )  
       call massdair(p,masse)  
       CALL exner_hyb(psi,p,pks,pk,pkf)  
   
 c    ....  Calcul de pls , pression au milieu des couches ,en Pascals  
       unskap=1./kappa  
       prefkap =  preff  ** kappa  
 c     PRINT *,' Pref kappa unskap  ',preff,kappa,unskap  
       DO l = 1, llm  
        DO j=1,jjp1  
         DO i =1, iip1  
         pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap  
         ENDDO  
        ENDDO  
        ENDDO  
   
   
 c -----------------------------------------------------------------  
 c   calcul des pressions pour les grilles u et v  
 c -----------------------------------------------------------------  
   
       do l=1,llm  
       do j=1,jjp1  
          do i=1,iip1  
             pext(i,j,l)=pls(i,j,l)*aire_2d(i,j)  
          enddo  
       enddo  
       enddo  
       call massbar(pext, pbarx, pbary )  
       do l=1,llm  
       do j=1,jjp1  
          do i=1,iip1  
             plunc(i,jjp1+1-j,l)=pbarx(i,j,l)/aireu_2d(i,j)  
             plsnc(i,jjp1+1-j,l)=pls(i,j,l)  
          enddo  
       enddo  
       enddo  
       do l=1,llm  
       do j=1,jjm  
          do i=1,iip1  
             plvnc(i,jjm+1-j,l)=pbary(i,j,l)/airev_2d(i,j)  
          enddo  
       enddo  
       enddo  
   
 c -----------------------------------------------------------------  
   
       if (guide_P) then  
       do j=1,jjp1  
          do i=1,iim  
             ps(i,j)=psnc(i,jjp1+1-j)  
          enddo  
          ps(iip1,j)=ps(1,j)  
       enddo  
       endif  
   
   
 c -----------------------------------------------------------------  
       call pres2lev(unc,zu,nlevnc,llm,pl,plunc,iip1,jjp1)  
       call pres2lev(vnc,zv,nlevnc,llm,pl,plvnc,iip1,jjm )  
       call pres2lev(tnc,zt,nlevnc,llm,pl,plsnc,iip1,jjp1)  
       call pres2lev(qnc,zq,nlevnc,llm,pl,plsnc,iip1,jjp1)  
   
 c     call dump2d(iip1,jjp1,ps,'PS    ')  
 c     call dump2d(iip1,jjp1,psu,'PS    ')  
 c     call dump2d(iip1,jjm,psv,'PS    ')  
 c  Inversion Nord/Sud  
       do l=1,llm  
          do j=1,jjp1  
             do i=1,iim  
                u(i,j,l)=zu(i,jjp1+1-j,l)  
                t(i,j,l)=zt(i,jjp1+1-j,l)  
                q(i,j,l)=zq(i,jjp1+1-j,l)  
             enddo  
             u(iip1,j,l)=u(1,j,l)  
             t(iip1,j,l)=t(1,j,l)  
             q(iip1,j,l)=q(1,j,l)  
          enddo  
       enddo  
   
       do l=1,llm  
          do j=1,jjm  
             do i=1,iim  
                v(i,j,l)=zv(i,jjm+1-j,l)  
             enddo  
             v(iip1,j,l)=v(1,j,l)  
          enddo  
       enddo  
   
       return  
       end  
   
 c===========================================================================  
       subroutine nat2gcm(u,v,t,rh,pk,ucov,vcov,teta,q)  
 c===========================================================================  
   
       use dimens_m  
       use paramet_m  
       use comconst  
       use comvert  
       use comgeom  
       use q_sat_m, only: q_sat  
       use guide_m  
       implicit none  
   
   
       real u(iip1,jjp1,llm),v(iip1,jjm,llm)  
       real t(iip1,jjp1,llm),pk(iip1,jjp1,llm),rh(iip1,jjp1,llm)  
       real ps(iip1,jjp1)  
   
       real ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm)  
       real teta(iip1,jjp1,llm),q(iip1,jjp1,llm)  
   
       real pres(iip1,jjp1,llm),qsat(iip1,jjp1,llm)  
   
       real unskap  
   
       integer i,j,l  
   
   
       print*,'Entree dans nat2gcm'  
 c    ucov(:,:,:)=0.  
 c    do l=1,llm  
 c       ucov(:,2:jjm,l)=u(:,2:jjm,l)*cu_2d(:,2:jjm)  
 c    enddo  
 c    ucov(iip1,:,:)=ucov(1,:,:)  
   
 c    teta(:,:,:)=t(:,:,:)*cpp/pk(:,:,:)  
 c    teta(iip1,:,:)=teta(1,:,:)  
       
 c   calcul de ucov et de la temperature potentielle  
       do l=1,llm  
          do j=1,jjp1  
             do i=1,iim  
                ucov(i,j,l)=u(i,j,l)*cu_2d(i,j)  
                teta(i,j,l)=t(i,j,l)*cpp/pk(i,j,l)  
             enddo  
             ucov(iip1,j,l)=ucov(1,j,l)  
             teta(iip1,j,l)=teta(1,j,l)  
          enddo  
          do i=1,iip1  
             ucov(i,1,l)=0.  
             ucov(i,jjp1,l)=0.  
             teta(i,1,l)=teta(1,1,l)  
             teta(i,jjp1,l)=teta(1,jjp1,l)  
          enddo  
       enddo  
   
 c   calcul de ucov  
       do l=1,llm  
          do j=1,jjm  
             do i=1,iim  
                vcov(i,j,l)=v(i,j,l)*cv_2d(i,j)  
             enddo  
             vcov(iip1,j,l)=vcov(1,j,l)  
          enddo  
       enddo  
   
 c     call dump2d(iip1,jjp1,teta,'TETA EN BAS   ')  
 c     call dump2d(iip1,jjp1,teta(1,1,llm),'TETA EN HAUT   ')  
   
 c  Humidite relative -> specifique  
 c  -------------------------------  
       if (1.eq.0) then  
 c   FINALEMENT ON GUIDE EN HUMIDITE RELATIVE  
       print*,'calcul de unskap'  
       unskap   = 1./ kappa  
       print*,'calcul de pres'  
       pres(:,:,:)=preff*(pk(:,:,:)/cpp)**unskap  
       print*,'calcul de qsat'  
       qsat = q_sat(t, pres)  
       print*,'calcul de q'  
 c   ATTENTION : humidites relatives en %  
       rh(:,:,:)=max(rh(:,:,:)*0.01,1.e-6)  
       q(:,:,:)=qsat(:,:,:)*rh(:,:,:)  
       print*,'calcul de q OK'  
   
       call dump2d(iip1,jjp1,pres,'PRESSION PREMIERE COUCHE   ')  
       call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1   ')  
       endif  
   
   
       return  
       end  
   
   
   
 c===========================================================================  
       subroutine correctbid(iim,nl,x)  
 c===========================================================================  
       integer iim,nl  
       real x(iim+1,nl)  
       integer i,l  
       real zz  
   
       do l=1,nl  
          do i=2,iim-1  
             if(abs(x(i,l)).gt.1.e10) then  
                zz=0.5*(x(i-1,l)+x(i+1,l))  
 c              print*,'correction ',i,l,x(i,l),zz  
                x(i,l)=zz  
             endif  
          enddo  
       enddo  
       return  
       end  
Legend:
Removed from v.3  
changed lines
  Added in v.134
  ViewVC Help
Powered by ViewVC 1.1.21