/[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 172 by guez, Wed Sep 30 15:59:14 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(psi, u, v, t, q)
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_nowrite
14  c -----------------------------------------------------------------      USE netcdf95, ONLY: nf95_get_var, nf95_inq_dimid, nf95_inq_varid, &
15        use dimens_m           nf95_inquire_dimension, nf95_open
16        use paramet_m      USE paramet_m, ONLY: iip1, jjp1
17        use comvert      use reanalyse2nat_m, only: reanalyse2nat
18        use comgeom  
19        use guide_m      real, intent(in):: psi(:, :) ! (iip1, jjp1)
20        IMPLICIT NONE      real, intent(out):: u(:, :, :) ! (iip1, jjp1, llm)
21        real, intent(out):: v(:, :, :) ! (iip1, jjm, llm)
22  c common      real, intent(out):: t(:, :, :), q(:, :, :) ! (iip1, jjp1, llm)
23  c ------  
24        ! Local:
25        include "netcdf.inc"      integer, save:: nlevnc
26        integer:: timestep = 0
27        real pk(iip1, jjp1, llm)
28  c arguments      integer, save:: ncidu, varidu, ncidv, varidv, ncidt, varidt, ncidQ, varidQ
29  c ---------      integer ncid, varid, dimid
30        integer nlevnc      real, allocatable, save:: unc(:, :, :) ! (iip1, jjp1, nlevnc)
31        integer timestep,mode,l      real, allocatable, save:: vnc(:, :, :) ! (iip1, jjm, nlevnc)
32        real, allocatable, save:: tnc(:, :, :), Qnc(:, :, :) ! (iip1, jjp1, nlevnc)
33        real psi(iip1,jjp1)      real, allocatable, save:: pl(:) ! (nlevnc)
34        real u(iip1,jjp1,llm),v(iip1,jjm,llm)      logical:: first = .true.
35        real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)      character(len = 8) name
36        real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)  
37        ! -----------------------------------------------------------------
38    
39  c local      ! Initialisation de la lecture des fichiers
40  c -----  
41        integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps      if (first) then
42        integer ncidpl         print *, 'Intitialisation de read reanalsye'
43        integer varidpl,ncidQ,varidQ  
44        save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps         ! Vent zonal
45        save ncidpl         if (guide_u) then
46        save varidpl,ncidQ,varidQ            call nf95_open('u.nc', nf90_nowrite, ncidu)
47              call nf95_inq_varid(ncidu, 'UWND', varidu)
48        real*4 unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)         endif
49        real*4 tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)  
50        real*4 Qnc(iip1,jjp1,nlevnc)         ! Vent meridien
51        real*4 pl(nlevnc)         if (guide_v) then
52              call nf95_open('v.nc', nf90_nowrite, ncidv)
53        integer start(4),count(4),status            call nf95_inq_varid(ncidv, 'VWND', varidv)
54           endif
55        real rcode  
56        logical first         ! Temperature
57        save first         if (guide_T) then
58              call nf95_open('T.nc', nf90_nowrite, ncidt)
59        data first/.true./            call nf95_inq_varid(ncidt, 'AIR', varidt)
60           endif
61    
62           ! Humidite
63  c -----------------------------------------------------------------         if (guide_Q) then
64  c   Initialisation de la lecture des fichiers            call nf95_open('hur.nc', nf90_nowrite, ncidQ)
65  c -----------------------------------------------------------------            call nf95_inq_varid(ncidQ, 'RH', varidQ)
66        if (first) then         endif
67             ncidpl=-99  
68             print*,'Intitialisation de read reanalsye'         ! Coordonn\'ee verticale :
69    
70  c Vent zonal         if (guide_u) then
71              if (guide_u) then            ncid = ncidu
72              ncidu=NCOPN('u.nc',NCNOWRIT,rcode)         else if (guide_v) then
73              varidu=NCVID(ncidu,'UWND',rcode)            ncid = ncidv
74              print*,'ncidu,varidu',ncidu,varidu         else if (guide_T) then
75              if (ncidpl.eq.-99) ncidpl=ncidu            ncid = ncidt
76              endif         else
77              ncid = ncidq
78  c Vent meridien         end if
79              if (guide_v) then  
80              ncidv=NCOPN('v.nc',NCNOWRIT,rcode)         name = merge('LEVEL   ', 'PRESSURE', ncep)
81              varidv=NCVID(ncidv,'VWND',rcode)         call nf95_inq_dimid(ncid, name, dimid)
82              print*,'ncidv,varidv',ncidv,varidv         call nf95_inquire_dimension(ncid, dimid, nclen = nlevnc)
83              if (ncidpl.eq.-99) ncidpl=ncidv         call nf95_inq_varid(ncid, name, varid)
84              endif         PRINT *, 'nlevnc = ', nlevnc
85           allocate(unc(iip1, jjp1, nlevnc), vnc(iip1, jjm, nlevnc))
86  c Temperature         allocate(tnc(iip1, jjp1, nlevnc), Qnc(iip1, jjp1, nlevnc), pl(nlevnc))
87              if (guide_T) then         call NF95_GET_VAR(ncid, varid, pl)
88              ncidt=NCOPN('T.nc',NCNOWRIT,rcode)         pl = 100. * pl ! passage en pascal
89              varidt=NCVID(ncidt,'AIR',rcode)         first = .false.
90              print*,'ncidt,varidt',ncidt,varidt      endif
91              if (ncidpl.eq.-99) ncidpl=ncidt  
92              endif      ! lecture des champs u, v, T
93    
94  c Humidite      timestep = timestep + 1
95              if (guide_Q) then      unc = 0.
96              ncidQ=NCOPN('hur.nc',NCNOWRIT,rcode)      vnc = 0.
97              varidQ=NCVID(ncidQ,'RH',rcode)      tnc = 0.
98              print*,'ncidQ,varidQ',ncidQ,varidQ      Qnc = 0.
99              if (ncidpl.eq.-99) ncidpl=ncidQ  
100              endif      ! Vent zonal
101        if (guide_u) then
102  c Pression de surface         call NF95_GET_VAR(ncidu, varidu, unc, start = (/1, 1, 1, timestep/))
103              if (guide_P) then         call correctbid(iim, jjp1 * nlevnc, unc)
104              ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)      endif
105              varidps=NCVID(ncidps,'SP',rcode)  
106              print*,'ncidps,varidps',ncidps,varidps      ! Temperature
107              endif      if (guide_T) then
108           call NF95_GET_VAR(ncidt, varidt, tnc, start = (/1, 1, 1, timestep/))
109  c Coordonnee verticale         call correctbid(iim, jjp1 * nlevnc, tnc)
110              if (ncep) then      endif
111                 print*,'Vous etes entrain de lire des donnees NCEP'  
112                 varidpl=NCVID(ncidpl,'LEVEL',rcode)      ! Humidite
113              else      if (guide_Q) then
114                 print*,'Vous etes entrain de lire des donnees ECMWF'         call NF95_GET_VAR(ncidQ, varidQ, Qnc, start = (/1, 1, 1, timestep/))
115                 varidpl=NCVID(ncidpl,'PRESSURE',rcode)         call correctbid(iim, jjp1 * nlevnc, Qnc)
116              endif      endif
117              print*,'ncidu,varidpl',ncidu,varidpl  
118        endif      ! Vent meridien
119        print*,'ok1'      if (guide_v) then
120           call NF95_GET_VAR(ncidv, varidv, vnc, start = (/1, 1, 1, timestep/))
121  c Niveaux de pression         call correctbid(iim, jjm * nlevnc, vnc)
122        print*,'WARNING!!! Il n y a pas de test de coherence'      endif
123        print*,'sur le nombre de niveaux verticaux dans le fichier nc'  
124        status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,pl)      call reanalyse2nat(nlevnc, psi, unc, vnc, tnc, Qnc, pl, u, v, t, Q, pk)
125  c  passage en pascal      call nat2gcm(u, v, t, Q, pk, u, v, t, Q)
126        pl(:)=100.*pl(:)  
127        if (first) then    end subroutine read_reanalyse
128         do l=1,nlevnc  
129            print*,'PL(',l,')=',pl(l)  end module read_reanalyse_m
        enddo  
       endif  
   
 c -----------------------------------------------------------------  
 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.172
  ViewVC Help
Powered by ViewVC 1.1.21