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/dyn3dpar/read_reanalyse.F

    r764 r1046  
    2626#include "guide.h"
    2727
    28 
    2928c arguments
    3029c ---------
    3130      integer nlevnc
    32       integer timestep,mode,l
     31      integer timestep,mode
    3332
    3433      real psi(iip1,jjp1)
    … …  
    3736      real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
    3837
    39 
    4038c local
    4139c -----
    4240      integer ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
    4341      integer ncidpl
    44       integer varidpl,ncidQ,varidQ
    45       save ncidpl,ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
     42      integer varidpl,ncidQ,varidQ,varidap,varidbp
     43      save ncidu,varidu,ncidv,varidv,ncidt,varidt,ncidps,varidps
     44      save ncidpl,varidap,varidbp
    4645      save varidpl,ncidQ,varidQ
    4746
    … …  
    4948      real*4 tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
    5049      real*4 Qnc(iip1,jjp1,nlevnc)
    51       real*4 pl(nlevnc),presnc(iip1,jjp1,nlevnc)
     50      real*4 plnc(iip1,jjp1,nlevnc)
     51      real*4 apnc(nlevnc),bpnc(nlevnc)
    5252
    5353      integer start(4),count(4),status
     54      integer i,j,l
    5455
    5556      real rcode
    … …  
    6061      data first/.true./
    6162
    62 
    63 
    6463c -----------------------------------------------------------------
    6564c   Initialisation de la lecture des fichiers
    … …  
    6867           ncidpl=-99
    6968           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
    7079
    7180c Vent zonal
    … …  
    8291               varidv=NCVID(ncidv,'VWND',rcode)
    8392               print*,'ncidv,varidv',ncidv,varidv
    84                if (ncidpl.eq.-99) ncidpl=ncidu
    85             endif
    86            
     93               if (ncidpl.eq.-99) ncidpl=ncidv
     94            endif
    8795
    8896c Temperature
    … …  
    9199               varidt=NCVID(ncidt,'AIR',rcode)
    92100               print*,'ncidt,varidt',ncidt,varidt
    93                if (ncidpl.eq.-99) ncidpl=ncidu
    94             endif
    95            
     101               if (ncidpl.eq.-99) ncidpl=ncidt
     102            endif
     103
    96104c Humidite
    97105            if (guide_Q) then
    … …  
    99107               varidQ=NCVID(ncidQ,'RH',rcode)
    100108               print*,'ncidQ,varidQ',ncidQ,varidQ
    101                if (ncidpl.eq.-99) ncidpl=ncidu
    102             endif
    103            
     109               if (ncidpl.eq.-99) ncidpl=ncidQ
     110            endif
     111
    104112c Pression de surface
    105             if (guide_P) then
     113            if ((guide_P).OR.(guide_modele)) then
    106114               ncidps=NCOPN('ps.nc',NCNOWRIT,rcode)
    107115               varidps=NCVID(ncidps,'SP',rcode)
    108116               print*,'ncidps,varidps',ncidps,varidps
    109117            endif
    110            
    111 c Coordonnee vertcale
    112             if (ncep) then
     118
     119c Coordonnee verticale
     120            if (.not.guide_modele) then
     121              if (ncep) then
    113122               print*,'Vous etes entrain de lire des donnees NCEP'
    114123               varidpl=NCVID(ncidpl,'LEVEL',rcode)
    … …  
    116125               print*,'Vous etes entrain de lire des donnees ECMWF'
    117126               varidpl=NCVID(ncidpl,'PRESSURE',rcode)
     127              endif
     128              print*,'ncidpl,varidpl',ncidpl,varidpl
    118129            endif
    119130            print*,'ncidu,varidpl',ncidu,varidpl
    120131
    121132      endif
    122 
    123133      print*,'ok1'
    124134
    … …  
    126136c   lecture des champs u, v, T, ps
    127137c -----------------------------------------------------------------
    128 
    129 c  niveaux de pression
    130 c  -------------------
    131 
    132       print*,'WARNING!!! Il n y a pas de test de coherence'
    133       print*,'sur le nombre de niveaux verticaux dans le fichier nc'
    134 #ifdef NC_DOUBLE
    135       status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,pl)
    136 #else
    137       status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,pl)
    138 #endif
    139      
    140 c  passage en pascal
    141       pl(:)=100.*pl(:)
    142 c  passage des presions sur une grille 3D
    143       do l=1,nlevnc
    144          presnc(:,:,l)=pl(l)
    145       enddo
    146       if(first) then
    147          do l=1,nlevnc
    148             print*,'PL(',l,')=',pl(l)
    149          enddo
    150       endif
    151 
    152138
    153139c  dimensions pour les champs scalaires et le vent zonal
    … …  
    163149      count(3)=nlevnc
    164150      count(4)=1
    165 
    166151
    167152c mise a zero des tableaux
    … …  
    171156       tnc(:,:,:)=0.
    172157       Qnc(:,:,:)=0.
     158       plnc(:,:,:)=0.
    173159
    174160c  Vent zonal
    … …  
    183169         status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unc)
    184170#endif
    185 c     call dump2d(iip1,jjp1,unc,'VENT NCEP   ')
    186 c     call dump2d(iip1,40,unc(1,1,nlevnc),'VENT NCEP   ')
    187          print*,'WARNING!!! Correction bidon pour palier a un '
    188          print*,'probleme dans la creation des fichiers nc'
    189          call correctbid(iim,jjp1*nlevnc,unc)
    190          call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ')
     171      print*,'WARNING!!! Correction bidon pour palier a un '
     172      print*,'probleme dans la creation des fichiers nc'
     173      call correctbid(iim,jjp1*nlevnc,unc)
     174      call dump2d(iip1,jjp1,unc,'UNC COUCHE 1 ')
    191175      endif
    192176
    … …  
    242226c  Pression de surface
    243227c  -------------------
    244 
    245       if (guide_P) then
     228      psnc(:,:)=0.
     229      if ((guide_P).OR.(guide_modele)) then
    246230#ifdef NC_DOUBLE
    247231         status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnc)
    … …  
    253237      endif
    254238
    255 
     239c Calcul de la pression aux differents niveaux
     240c --------------------------------------------
     241      if (guide_modele) then
     242#ifdef NC_DOUBLE
     243      status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
     244      status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
     245#else
     246      status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
     247      status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
     248#endif
     249      else
     250#ifdef NC_DOUBLE
     251      status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
     252#else
     253      status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
     254#endif
     255c conversion en Pascals
     256      apnc=apnc*100.
     257      bpnc(:)=0.
     258      endif
     259      do i=1,iip1
     260        do j=1,jjp1
     261          do l=1,nlevnc
     262            plnc(i,j,l)=apnc(l)+bpnc(l)*psnc(i,j)
     263          enddo
     264        enddo
     265      enddo
     266      if (first) then
     267        print*,'Verification ordre niveaux verticaux'
     268        print*,'LMDZ :'
     269        do l=1,llm
     270           print*,'PL(',l,')=',(ap(l)+ap(l+1))/2.
     271     $                    +psi(1,jjp1)*(bp(l)+bp(l+1))/2.
     272        enddo
     273        print*,'Fichiers guidage'
     274        do l=1,nlevnc
     275           print*,'PL(',l,')=',plnc(1,1,l)
     276        enddo
     277        if (guide_u) then
     278          do l=1,nlevnc
     279            print*,'U(',l,')=',unc(1,1,l)
     280          enddo
     281        endif
     282        if (guide_T) then
     283          do l=1,nlevnc
     284            print*,'T(',l,')=',tnc(1,1,l)
     285          enddo
     286        endif
     287        print *,'inversion de l''ordre: ok_invertp=',ok_invertp
     288      endif
    256289
    257290c -----------------------------------------------------------------
    258291c  Interpollation verticale sur les niveaux modele
    259292c -----------------------------------------------------------------
    260       call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,pl,u,v,t,Q
     293      call reanalyse2nat(nlevnc,psi,unc,vnc,tnc,Qnc,psnc,plnc,u,v,t,Q
    261294     s    ,ps,masse,pk)
    262295
    … …  
    288321c===========================================================================
    289322      subroutine reanalyse2nat(nlevnc,psi
    290      s   ,unc,vnc,tnc,qnc,psnc,pl,u,v,t,q
     323     s   ,unc,vnc,tnc,qnc,psnc,plnc,u,v,t,q
    291324     s   ,ps,masse,pk)
    292325c===========================================================================
    … …  
    311344      real t(iip1,jjp1,llm),ps(iip1,jjp1),q(iip1,jjp1,llm)
    312345
    313       real pl(nlevnc)
     346      real plnc(iip1,jjp1,nlevnc)
    314347      real unc(iip1,jjp1,nlevnc),vnc(iip1,jjm,nlevnc)
    315348      real tnc(iip1,jjp1,nlevnc),psnc(iip1,jjp1)
    … …  
    396429
    397430c -----------------------------------------------------------------
    398       call pres2lev(unc,zu,nlevnc,llm,pl,plunc,iip1,jjp1)
    399       call pres2lev(vnc,zv,nlevnc,llm,pl,plvnc,iip1,jjm )
    400       call pres2lev(tnc,zt,nlevnc,llm,pl,plsnc,iip1,jjp1)
    401       call pres2lev(qnc,zq,nlevnc,llm,pl,plsnc,iip1,jjp1)
     431      call pres2lev(unc,zu,nlevnc,llm,plnc,plunc,iip1,jjp1,ok_invertp)
     432      call pres2lev(vnc,zv,nlevnc,llm,plnc,plvnc,iip1,jjm,ok_invertp )
     433      call pres2lev(tnc,zt,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp)
     434      call pres2lev(qnc,zq,nlevnc,llm,plnc,plsnc,iip1,jjp1,ok_invertp)
    402435
    403436c     call dump2d(iip1,jjp1,ps,'PS    ')
    … …  
    500533c  Humidite relative -> specifique
    501534c  -------------------------------
    502       if (1.eq.0) then
     535      if (guide_hr) then
    503536c   FINALEMENT ON GUIDE EN HUMIDITE RELATIVE
    504537      print*,'calcul de unskap'
    … …  
    513546      q(:,:,:)=qsat(:,:,:)*rh(:,:,:)
    514547      print*,'calcul de q OK'
    515 
     548      call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1   ')
    516549      call dump2d(iip1,jjp1,pres,'PRESSION PREMIERE COUCHE   ')
     550      else
     551      q(:,:,:)=rh(:,:,:)
     552      print*,'calcul de q OK'
    517553      call dump2d(iip1,jjp1,q,'HUMIDITE SPECIFIQUE COUCHE 1   ')
    518       endif
    519 
     554      endif
    520555
    521556      return
Note: See TracChangeset for help on using the changeset viewer.