/[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 20 by guez, Wed Oct 15 16:19:57 2008 UTC trunk/libf/dyn3d/Read_reanalyse/read_reanalyse.f90 revision 37 by guez, Tue Dec 21 15:45:48 2010 UTC
# Line 1  Line 1 
1          subroutine read_reanalyse(timestep,psi &
2             ,u,v,t,q,masse,ps,mode,nlevnc)
3    
4  !  !
5  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/read_reanalyse.F,v 1.3 2005/04/15 12:31:21 lmdzadmin Exp $  ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/read_reanalyse.F,v 1.3 2005/04/15 12:31:21 lmdzadmin Exp $
6  !  !
7  c  !
8  c  !
9        subroutine read_reanalyse(timestep,psi  !   mode=0 variables naturelles
10       s   ,u,v,t,q,masse,ps,mode,nlevnc)  !   mode=1 variabels GCM
11    
12  c   mode=0 variables naturelles  ! -----------------------------------------------------------------
13  c   mode=1 variabels GCM  !   Declarations
14    ! -----------------------------------------------------------------
 c -----------------------------------------------------------------  
 c   Declarations  
 c -----------------------------------------------------------------  
15        use dimens_m        use dimens_m
16        use paramet_m        use paramet_m
17        use comvert        use comvert
# Line 21  c -------------------------------------- Line 21  c --------------------------------------
21    
22        IMPLICIT NONE        IMPLICIT NONE
23    
24  c common  ! common
25  c ------  ! ------
26    
27        include "netcdf.inc"        include "netcdf.inc"
28    
29    
30  c arguments  ! arguments
31  c ---------  ! ---------
32        integer nlevnc        integer nlevnc
33        integer timestep,mode,l        integer timestep,mode,l
34    
# Line 38  c --------- Line 38  c ---------
38        real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)        real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
39    
40    
41  c local  ! local
42  c -----  ! -----
43        integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps        integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
44        integer ncidpl        integer ncidpl
45        integer varidpl,ncidQ,varidQ        integer varidpl,ncidQ,varidQ
# Line 47  c ----- Line 47  c -----
47        save ncidpl        save ncidpl
48        save varidpl,ncidQ,varidQ        save varidpl,ncidQ,varidQ
49    
50        real*4 unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)        real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
51        real*4 tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)        real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
52        real*4 Qnc(iip1,jjp1,nlevnc)        real Qnc(iip1,jjp1,nlevnc)
53        real*4 pl(nlevnc)        real pl(nlevnc)
54    
55        integer start(4),count(4),status        integer start(4),count(4),status
56    
# Line 62  c ----- Line 62  c -----
62    
63    
64    
65  c -----------------------------------------------------------------  ! -----------------------------------------------------------------
66  c   Initialisation de la lecture des fichiers  !   Initialisation de la lecture des fichiers
67  c -----------------------------------------------------------------  ! -----------------------------------------------------------------
68        if (first) then        if (first) then
69             ncidpl=-99             ncidpl=-99
70             print*,'Intitialisation de read reanalsye'             print*,'Intitialisation de read reanalsye'
71    
72  c Vent zonal  ! Vent zonal
73              if (guide_u) then              if (guide_u) then
74              rcode=nf90_open('u.nc',nf90_nowrite,ncidu)              rcode=nf90_open('u.nc',nf90_nowrite,ncidu)
75              rcode = nf90_inq_varid(ncidu, 'UWND', varidu)              rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
# Line 77  c Vent zonal Line 77  c Vent zonal
77              if (ncidpl.eq.-99) ncidpl=ncidu              if (ncidpl.eq.-99) ncidpl=ncidu
78              endif              endif
79    
80  c Vent meridien  ! Vent meridien
81              if (guide_v) then              if (guide_v) then
82              rcode=nf90_open('v.nc',nf90_nowrite,ncidv)              rcode=nf90_open('v.nc',nf90_nowrite,ncidv)
83              rcode = nf90_inq_varid(ncidv, 'VWND', varidv)              rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
# Line 85  c Vent meridien Line 85  c Vent meridien
85              if (ncidpl.eq.-99) ncidpl=ncidv              if (ncidpl.eq.-99) ncidpl=ncidv
86              endif              endif
87    
88  c Temperature  ! Temperature
89              if (guide_T) then              if (guide_T) then
90              rcode=nf90_open('T.nc',nf90_nowrite,ncidt)              rcode=nf90_open('T.nc',nf90_nowrite,ncidt)
91              rcode = nf90_inq_varid(ncidt, 'AIR', varidt)              rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
# Line 93  c Temperature Line 93  c Temperature
93              if (ncidpl.eq.-99) ncidpl=ncidt              if (ncidpl.eq.-99) ncidpl=ncidt
94              endif              endif
95    
96  c Humidite  ! Humidite
97              if (guide_Q) then              if (guide_Q) then
98              rcode=nf90_open('hur.nc',nf90_nowrite,ncidQ)              rcode=nf90_open('hur.nc',nf90_nowrite,ncidQ)
99              rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)              rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
# Line 101  c Humidite Line 101  c Humidite
101              if (ncidpl.eq.-99) ncidpl=ncidQ              if (ncidpl.eq.-99) ncidpl=ncidQ
102              endif              endif
103    
104  c Pression de surface  ! Pression de surface
105              if (guide_P) then              if (guide_P) then
106              rcode=nf90_open('ps.nc',nf90_nowrite,ncidps)              rcode=nf90_open('ps.nc',nf90_nowrite,ncidps)
107              rcode = nf90_inq_varid(ncidps, 'SP', varidps)              rcode = nf90_inq_varid(ncidps, 'SP', varidps)
108              print*,'ncidps,varidps',ncidps,varidps              print*,'ncidps,varidps',ncidps,varidps
109              endif              endif
110    
111  c Coordonnee verticale  ! Coordonnee verticale
112              if (ncep) then              if (ncep) then
113                 print*,'Vous etes entrain de lire des donnees NCEP'                 print*,'Vous etes entrain de lire des donnees NCEP'
114                 rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)                 rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
# Line 120  c Coordonnee verticale Line 120  c Coordonnee verticale
120        endif        endif
121        print*,'ok1'        print*,'ok1'
122    
123  c Niveaux de pression  ! Niveaux de pression
124        print*,'WARNING!!! Il n y a pas de test de coherence'        print*,'WARNING!!! Il n y a pas de test de coherence'
125        print*,'sur le nombre de niveaux verticaux dans le fichier nc'        print*,'sur le nombre de niveaux verticaux dans le fichier nc'
126        status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,pl)        status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,pl)
127  c  passage en pascal  !  passage en pascal
128        pl(:)=100.*pl(:)        pl(:)=100.*pl(:)
129        if (first) then        if (first) then
130         do l=1,nlevnc         do l=1,nlevnc
# Line 132  c  passage en pascal Line 132  c  passage en pascal
132         enddo         enddo
133        endif        endif
134    
135  c -----------------------------------------------------------------  ! -----------------------------------------------------------------
136  c   lecture des champs u, v, T, ps  !   lecture des champs u, v, T, ps
137  c -----------------------------------------------------------------  ! -----------------------------------------------------------------
138    
139  c  dimensions pour les champs scalaires et le vent zonal  !  dimensions pour les champs scalaires et le vent zonal
140  c  -----------------------------------------------------  !  -----------------------------------------------------
141    
142        start(1)=1        start(1)=1
143        start(2)=1        start(2)=1
# Line 149  c  ------------------------------------- Line 149  c  -------------------------------------
149        count(3)=nlevnc        count(3)=nlevnc
150        count(4)=1        count(4)=1
151    
152  c mise a zero des tableaux  ! mise a zero des tableaux
153  c ------------------------  ! ------------------------
154         unc(:,:,:)=0.         unc(:,:,:)=0.
155         vnc(:,:,:)=0.         vnc(:,:,:)=0.
156         tnc(:,:,:)=0.         tnc(:,:,:)=0.
157         Qnc(:,:,:)=0.         Qnc(:,:,:)=0.
158    
159  c  Vent zonal  !  Vent zonal
160  c  ----------  !  ----------
161    
162        if (guide_u) then        if (guide_u) then
163        print*,'avant la lecture de UNCEP nd de niv:',nlevnc        print*,'avant la lecture de UNCEP nd de niv:',nlevnc
164        status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)        status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)
165  c     call dump2d(iip1,jjp1,unc,'VENT NCEP   ')  !     call dump2d(iip1,jjp1,unc,'VENT NCEP   ')
166  c     call dump2d(iip1,40,unc(1,1,nlevnc),'VENT NCEP   ')  !     call dump2d(iip1,40,unc(1,1,nlevnc),'VENT NCEP   ')
167        print*,'WARNING!!! Correction bidon pour palier a un '        print*,'WARNING!!! Correction bidon pour palier a un '
168        print*,'probleme dans la creation des fichiers nc'        print*,'probleme dans la creation des fichiers nc'
169        call correctbid(iim,jjp1*nlevnc,unc)        call correctbid(iim,jjp1*nlevnc,unc)
170        call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ')        call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ')
171        endif        endif
172    
173  c  Temperature  !  Temperature
174  c  -----------  !  -----------
175    
176        print*,'ncidt=',ncidt,'varidt=',varidt,'start=',start        print*,'ncidt=',ncidt,'varidt=',varidt,'start=',start
177        print*,'count=',count        print*,'count=',count
# Line 182  c  ----------- Line 182  c  -----------
182        call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 BBB ')        call dump2d(iip1,jjp1,tnc,'TNC COUCHE 1 BBB ')
183        endif        endif
184    
185  c  Humidite  !  Humidite
186  c  --------  !  --------
187    
188        if (guide_Q) then        if (guide_Q) then
189        status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,Qnc)        status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,Qnc)
# Line 192  c  -------- Line 192  c  --------
192        endif        endif
193    
194        count(2)=jjm        count(2)=jjm
195  c  Vent meridien  !  Vent meridien
196  c  -------------  !  -------------
197    
198        if (guide_v) then        if (guide_v) then
199        status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc)        status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnc)
# Line 207  c  ------------- Line 207  c  -------------
207        count(3)=1        count(3)=1
208        count(4)=0        count(4)=0
209    
210  c  Pression de surface  !  Pression de surface
211  c  -------------------  !  -------------------
212    
213        if (guide_P) then        if (guide_P) then
214        status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)        status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnc)
# Line 218  c  ------------------- Line 218  c  -------------------
218    
219    
220    
221  c -----------------------------------------------------------------  ! -----------------------------------------------------------------
222  c  Interpollation verticale sur les niveaux modele  !  Interpollation verticale sur les niveaux modele
223  c -----------------------------------------------------------------  ! -----------------------------------------------------------------
224        call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,pl,u,v,t,Q        call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,pl,u,v,t,Q &
225       s    ,ps,masse,pk)            ,ps,masse,pk)
226    
227        call dump2d(iip1,jjm,v,'V COUCHE APRES ')        call dump2d(iip1,jjm,v,'V COUCHE APRES ')
228    
229    
230  c -----------------------------------------------------------------  ! -----------------------------------------------------------------
231  c  Passage aux variables du modele (vents covariants, temperature  !  Passage aux variables du modele (vents covariants, temperature
232  c  potentielle et humidite specifique)  !  potentielle et humidite specifique)
233  c -----------------------------------------------------------------  ! -----------------------------------------------------------------
234        call nat2gcm(u,v,t,Q,pk,u,v,t,Q)        call nat2gcm(u,v,t,Q,pk,u,v,t,Q)
235        print*,'TIMESTEP ',timestep        print*,'TIMESTEP ',timestep
236        if(mode.ne.1) stop'mode pas egal 0'        if(mode.ne.1) stop'mode pas egal 0'
237  c     call dump2d(iip1,jjm,v,'VCOV COUCHE 1 ')  !     call dump2d(iip1,jjm,v,'VCOV COUCHE 1 ')
238    
239  c   Lignes introduites a une epoque pour un probleme oublie...  !   Lignes introduites a une epoque pour un probleme oublie...
240  c     do l=1,llm  !     do l=1,llm
241  c        do i=1,iip1  !        do i=1,iip1
242  c           v(i,1,l)=0.  !           v(i,1,l)=0.
243  c           v(i,jjm,l)=0.  !           v(i,jjm,l)=0.
244  c        enddo  !        enddo
245  c     enddo  !     enddo
246        first=.false.        first=.false.
247    
248        return        return
249        end        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.20  
changed lines
  Added in v.37
  ViewVC Help
Powered by ViewVC 1.1.21