Ignore:
Timestamp:
Nov 6, 2008, 5:39:24 PM (16 years ago)
Author:
lmdzadmin
Message:

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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/dyn3d/read_reanalyse.F

    r618 r1046  
    2929c ---------
    3030      integer nlevnc
    31       integer timestep,mode,l
     31      integer timestep,mode
    3232
    3333      real psi(iip1,jjp1)
    … …  
    4141      integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
    4242      integer ncidpl
    43       integer varidpl,ncidQ,varidQ
     43      integer varidpl,ncidQ,varidQ,varidap,varidbp
    4444      save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
    45       save ncidpl
     45      save ncidpl,varidap,varidbp
    4646      save varidpl,ncidQ,varidQ
    4747
    … …  
    4949      real*4 tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
    5050      real*4 Qnc(iip1,jjp1,nlevnc)
    51       real*4 pl(nlevnc)
     51      real*4 plnc(iip1,jjp1,nlevnc)
     52      real*4 apnc(nlevnc),bpnc(nlevnc)
    5253
    5354      integer start(4),count(4),status
     55      integer i,j,l
    5456
    5557      real rcode
    … …  
    5961      data first/.true./
    6062
    61 
    62 
    6363c -----------------------------------------------------------------
    6464c   Initialisation de la lecture des fichiers
    … …  
    6767           ncidpl=-99
    6868           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
    6979
    7080c Vent zonal
    … …  
    101111
    102112c Pression de surface
    103             if (guide_P) then
     113            if ((guide_P).OR.(guide_modele)) then
    104114            ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
    105115            varidps=NCVID(ncidps,'SP',rcode)
    … …  
    108118
    109119c Coordonnee verticale
    110             if (ncep) then
     120            if (.not.guide_modele) then
     121              if (ncep) then
    111122               print*,'Vous etes entrain de lire des donnees NCEP'
    112123               varidpl=NCVID(ncidpl,'LEVEL',rcode)
    113             else
     124              else
    114125               print*,'Vous etes entrain de lire des donnees ECMWF'
    115126               varidpl=NCVID(ncidpl,'PRESSURE',rcode)
    116             endif
    117             print*,'ncidu,varidpl',ncidu,varidpl
     127              endif
     128              print*,'ncidpl,varidpl',ncidpl,varidpl
     129            endif
     130! endif (first)
    118131      endif
    119132      print*,'ok1'
    120133
    121 c Niveaux de pression
    122       print*,'WARNING!!! Il n y a pas de test de coherence'
    123       print*,'sur le nombre de niveaux verticaux dans le fichier nc'
    124 #ifdef NC_DOUBLE
    125       status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,pl)
    126 #else
    127       status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,pl)
    128 #endif
    129 c  passage en pascal
    130       pl(:)=100.*pl(:)
    131       if (first) then
    132        do l=1,nlevnc
    133           print*,'PL(',l,')=',pl(l)
    134        enddo
    135       endif
    136134
    137135c -----------------------------------------------------------------
    … …  
    158156       tnc(:,:,:)=0.
    159157       Qnc(:,:,:)=0.
     158       plnc(:,:,:)=0.
    160159
    161160c  Vent zonal
    … …  
    169168      status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)
    170169#endif
    171 c     call dump2d(iip1,jjp1,unc,'VENT NCEP   ')
    172 c     call dump2d(iip1,40,unc(1,1,nlevnc),'VENT NCEP   ')
    173170      print*,'WARNING!!! Correction bidon pour palier a un '
    174171      print*,'probleme dans la creation des fichiers nc'
    … …  
    228225c  Pression de surface
    229226c  -------------------
    230 
    231       if (guide_P) then
     227      psnc(:,:)=0.
     228      if ((guide_P).OR.(guide_modele)) then
    232229#ifdef NC_DOUBLE
    233230      status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnc)
    … …  
    239236      endif
    240237
    241 
     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
    242288
    243289c -----------------------------------------------------------------
    244290c  Interpollation verticale sur les niveaux modele
    245291c -----------------------------------------------------------------
    246       call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,pl,u,v,t,Q
     292      call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,plnc,u,v,t,Q
    247293     s    ,ps,masse,pk)
    248294
    … …  
    274320c===========================================================================
    275321      subroutine reanalyse2nat(nlevnc,psi
    276      s   ,unc,vnc,tnc,qnc,psnc,pl,u,v,t,q
     322     s   ,unc,vnc,tnc,qnc,psnc,plnc,u,v,t,q
    277323     s   ,ps,masse,pk)
    278324c===========================================================================
    … …  
    297343      real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)
    298344
    299       real pl(nlevnc)
     345      real plnc(iip1,jjp1,nlevnc)
    300346      real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
    301347      real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
    … …  
    382428
    383429c -----------------------------------------------------------------
    384       call pres2lev(unc,zu,nlevnc,llm,pl,plunc,iip1,jjp1)
    385       call pres2lev(vnc,zv,nlevnc,llm,pl,plvnc,iip1,jjm )
    386       call pres2lev(tnc,zt,nlevnc,llm,pl,plsnc,iip1,jjp1)
    387       call pres2lev(qnc,zq,nlevnc,llm,pl,plsnc,iip1,jjp1)
     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)
    388434
    389435c     call dump2d(iip1,jjp1,ps,'PS    ')
    … …  
    486532c  Humidite relative -> specifique
    487533c  -------------------------------
    488       if (1.eq.0) then
     534      if (guide_hr) then
    489535c   FINALEMENT ON GUIDE EN HUMIDITE RELATIVE
    490536      print*,'calcul de unskap'
    … …  
    499545      q(:,:,:)=qsat(:,:,:)*rh(:,:,:)
    500546      print*,'calcul de q OK'
    501 
     547      call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1   ')
    502548      call dump2d(iip1,jjp1,pres,'PRESSION PREMIERE COUCHE   ')
     549      else
     550      q(:,:,:)=rh(:,:,:)
     551      print*,'calcul de q OK'
    503552      call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1   ')
    504       endif
    505 
     553      endif
    506554
    507555      return
Note: See TracChangeset for help on using the changeset viewer.