/[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/dyn3d/Read_reanalyse/read_reanalyse.f revision 102 by guez, Tue Jul 15 13:43:24 2014 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 dump2d_m, only: dump2d
14  c -----------------------------------------------------------------      USE netcdf, ONLY: nf90_get_var, nf90_inq_varid, nf90_nowrite, nf90_open
15        use dimens_m      USE paramet_m, ONLY: iip1, jjp1
16        use paramet_m      use reanalyse2nat_m, only: reanalyse2nat
17        use comvert  
18        use comgeom      integer timestep
19        use guide_m      real, intent(in):: psi(iip1, jjp1)
20        IMPLICIT NONE      real u(iip1, jjp1, llm), v(iip1, jjm, llm)
21        real t(iip1, jjp1, llm), q(iip1, jjp1, llm)
22  c common      real masse(iip1, jjp1, llm)
23  c ------      integer nlevnc
24    
25        include "netcdf.inc"      ! Local:
26    
27        integer l
28  c arguments      real pk(iip1, jjp1, llm)
29  c ---------      integer, save:: ncidu, varidu, ncidv, varidv, ncidt, varidt
30        integer nlevnc      integer, save:: ncidpl
31        integer timestep,mode,l      integer, save:: varidpl, ncidQ, varidQ
32        real unc(iip1, jjp1, nlevnc), vnc(iip1, jjm, nlevnc)
33        real psi(iip1,jjp1)      real tnc(iip1, jjp1, nlevnc)
34        real u(iip1,jjp1,llm),v(iip1,jjm,llm)      real Qnc(iip1, jjp1, nlevnc)
35        real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)      real pl(nlevnc)
36        real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)      integer start(4), count(4), status
37        real rcode
38        logical:: first = .true.
39  c local  
40  c -----      ! -----------------------------------------------------------------
41        integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps  
42        integer ncidpl      !   Initialisation de la lecture des fichiers
43        integer varidpl,ncidQ,varidQ  
44        save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps      if (first) then
45        save ncidpl         ncidpl=-99
46        save varidpl,ncidQ,varidQ         print *, 'Intitialisation de read reanalsye'
47    
48        real*4 unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)         ! Vent zonal
49        real*4 tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)         if (guide_u) then
50        real*4 Qnc(iip1,jjp1,nlevnc)            rcode=nf90_open('u.nc', nf90_nowrite, ncidu)
51        real*4 pl(nlevnc)            rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
52              print *, 'ncidu, varidu', ncidu, varidu
53        integer start(4),count(4),status            if (ncidpl.eq.-99) ncidpl=ncidu
54           endif
55        real rcode  
56        logical first         ! Vent meridien
57        save first         if (guide_v) then
58              rcode=nf90_open('v.nc', nf90_nowrite, ncidv)
59        data first/.true./            rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
60              print *, 'ncidv, varidv', ncidv, varidv
61              if (ncidpl.eq.-99) ncidpl=ncidv
62           endif
63  c -----------------------------------------------------------------  
64  c   Initialisation de la lecture des fichiers         ! Temperature
65  c -----------------------------------------------------------------         if (guide_T) then
66        if (first) then            rcode=nf90_open('T.nc', nf90_nowrite, ncidt)
67             ncidpl=-99            rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
68             print*,'Intitialisation de read reanalsye'            print *, 'ncidt, varidt', ncidt, varidt
69              if (ncidpl.eq.-99) ncidpl=ncidt
70  c Vent zonal         endif
71              if (guide_u) then  
72              ncidu=NCOPN('u.nc',NCNOWRIT,rcode)         ! Humidite
73              varidu=NCVID(ncidu,'UWND',rcode)         if (guide_Q) then
74              print*,'ncidu,varidu',ncidu,varidu            rcode=nf90_open('hur.nc', nf90_nowrite, ncidQ)
75              if (ncidpl.eq.-99) ncidpl=ncidu            rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
76              endif            print *, 'ncidQ, varidQ', ncidQ, varidQ
77              if (ncidpl.eq.-99) ncidpl=ncidQ
78  c Vent meridien         endif
79              if (guide_v) then  
80              ncidv=NCOPN('v.nc',NCNOWRIT,rcode)         ! Coordonnee verticale
81              varidv=NCVID(ncidv,'VWND',rcode)         if (ncep) then
82              print*,'ncidv,varidv',ncidv,varidv            print *, 'Vous etes entrain de lire des donnees NCEP'
83              if (ncidpl.eq.-99) ncidpl=ncidv            rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
84              endif         else
85              print *, 'Vous etes entrain de lire des donnees ECMWF'
86  c Temperature            rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
87              if (guide_T) then         endif
88              ncidt=NCOPN('T.nc',NCNOWRIT,rcode)         print *, 'ncidu, varidpl', ncidu, varidpl
89              varidt=NCVID(ncidt,'AIR',rcode)      endif
90              print*,'ncidt,varidt',ncidt,varidt      print *, 'ok1'
91              if (ncidpl.eq.-99) ncidpl=ncidt  
92              endif      ! Niveaux de pression
93        print *, 'WARNING!!! Il n y a pas de test de coherence'
94  c Humidite      print *, 'sur le nombre de niveaux verticaux dans le fichier nc'
95              if (guide_Q) then      status=NF90_GET_VAR(ncidpl, varidpl, pl)
96              ncidQ=NCOPN('hur.nc',NCNOWRIT,rcode)      !  passage en pascal
97              varidQ=NCVID(ncidQ,'RH',rcode)      pl(:)=100.*pl(:)
98              print*,'ncidQ,varidQ',ncidQ,varidQ      if (first) then
99              if (ncidpl.eq.-99) ncidpl=ncidQ         do l=1, nlevnc
100              endif            print *, 'PL(', l, ')=', pl(l)
   
 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)  
101         enddo         enddo
102        endif      endif
103    
104        !   lecture des champs u, v, T
105    
106        !  dimensions pour les champs scalaires et le vent zonal
107    
108        start(1)=1
109        start(2)=1
110        start(3)=1
111        start(4)=timestep
112    
113        count(1)=iip1
114        count(2)=jjp1
115        count(3)=nlevnc
116        count(4)=1
117    
118        ! mise a zero des tableaux
119    
120        unc(:, :, :)=0.
121        vnc(:, :, :)=0.
122        tnc(:, :, :)=0.
123        Qnc(:, :, :)=0.
124    
125        !  Vent zonal
126    
127        if (guide_u) then
128           print *, 'avant la lecture de UNCEP nd de niv:', nlevnc
129           status=NF90_GET_VAR(ncidu, varidu, unc, start, count)
130           print *, 'WARNING!!! Correction bidon pour palier a un '
131           print *, 'probleme dans la creation des fichiers nc'
132           call correctbid(iim, jjp1*nlevnc, unc)
133           call dump2d(iip1, jjp1, unc, 'UNC COUCHE 1 ')
134        endif
135    
136        !  Temperature
137    
138        print *, 'ncidt=', ncidt, 'varidt=', varidt, 'start=', start
139        print *, 'count=', count
140        if (guide_T) then
141           status=NF90_GET_VAR(ncidt, varidt, tnc, start, count)
142           call dump2d(iip1, jjp1, tnc, 'TNC COUCHE 1 AAA ')
143           call correctbid(iim, jjp1*nlevnc, tnc)
144           call dump2d(iip1, jjp1, tnc, 'TNC COUCHE 1 BBB ')
145        endif
146    
147        !  Humidite
148    
149        if (guide_Q) then
150           status=NF90_GET_VAR(ncidQ, varidQ, Qnc, start, count)
151           call correctbid(iim, jjp1*nlevnc, Qnc)
152           call dump2d(iip1, jjp1, Qnc, 'QNC COUCHE 1 ')
153        endif
154    
155        count(2)=jjm
156        !  Vent meridien
157    
158        if (guide_v) then
159           status=NF90_GET_VAR(ncidv, varidv, vnc, start, count)
160           call correctbid(iim, jjm*nlevnc, vnc)
161           call dump2d(iip1, jjm, vnc, 'VNC COUCHE 1 ')
162        endif
163    
164        start(3)=timestep
165        start(4)=0
166        count(2)=jjp1
167        count(3)=1
168        count(4)=0
169    
170        !  Interpolation verticale sur les niveaux modele
171    
172        call reanalyse2nat(nlevnc, psi, unc, vnc, tnc, Qnc, pl, u, v, t, Q, &
173             masse, pk)
174    
175        call dump2d(iip1, jjm, v, 'V COUCHE APRES ')
176    
177        !  Passage aux variables du modele (vents covariants, temperature
178        !  potentielle et humidite specifique)
179    
180        call nat2gcm(u, v, t, Q, pk, u, v, t, Q)
181        print *, 'TIMESTEP ', timestep
182        first=.false.
183    
184      end subroutine read_reanalyse
185    
186  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.102

  ViewVC Help
Powered by ViewVC 1.1.21