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

Location:
LMDZ4/trunk/libf/dyn3d
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r999 r1046  
    682682      ok_strato=.FALSE.
    683683      CALL getin('ok_strato',ok_strato)
     684
     685!Config  Key  = ok_gradsfile
     686!Config  Desc = activation des sorties grads du guidage
     687!Config  Def  = n
     688!Config  Help = active les sorties grads du guidage
     689
     690       ok_gradsfile = .FALSE.
     691       CALL getin('ok_gradsfile',ok_gradsfile)
    684692
    685693      write(lunout,*)' #########################################'
    … …  
    717725      write(lunout,*)' config_inca = ', config_inca
    718726      write(lunout,*)' ok_strato = ', ok_strato
     727      write(lunout,*)' ok_gradsfile = ', ok_gradsfile
    719728c
    720729      RETURN

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

    r524 r1046  
    1919      call getpar('online',1,online,'Index de controle du guide')
    2020      CALL getpar('ncep',.false.,ncep,'Coordonnee vert NCEP ou ECMWF')
     21      CALL getpar('guide_modele',.false.,guide_modele,
     22     $            'guidage sur niveaux modele')
    2123      CALL getpar('ini_anal',.false.,ini_anal,'Initial = analyse')
     24      CALL getpar('ok_invertp',.true.,ok_invertp,'niveaux p inverses')
    2225
    2326      CALL getpar('guide_u',.true.,guide_u,'guidage de u')
    … …  
    2629      CALL getpar('guide_P',.true.,guide_P,'guidage de P')
    2730      CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q')
     31      CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R')
     32      CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
    2833
    2934c   Constantes de rappel. Unite : fraction de jour
    … …  
    3843      CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
    3944      CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
     45      CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
    4046
    4147c   Latitude min et max pour le rappel.

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

    r1042 r1046  
    7878      real alpha_T(ip1jmp1),alpha_P(ip1jmp1)
    7979      real alpha_u(ip1jmp1),alpha_v(ip1jm)
     80      real alpha_pcor(llm)
    8081      real dday_step,toto,reste,itau_test
    8182      INTEGER step_rea,count_no_rea
    … …  
    111112
    112113      save alpha_T,alpha_q,alpha_u,alpha_v,alpha_P,itau_test
     114      save alpha_pcor
    113115      save step_rea,count_no_rea
    114116
    115 #define gradsfile
    116 #undef gradsfile
    117 
    118 #ifdef gradsfile
    119117      character*10 file
    120118      integer igrads
    … …  
    122120      save igrads,dtgrads
    123121      data igrads,dtgrads/2,100./
    124 #endif
    125122C-----------------------------------------------------------------------
    126123c calcul de l'humidite saturante
    … …  
    160157
    161158
    162 #ifdef gradsfile
     159        if (ok_gradsfile) then
    163160         file='guide'
    164161         call inigrads(igrads,iip1
    … …  
    166163     s  ,llm,presnivs,1.
    167164     s  ,dtgrads,file,'dyn_zon ')
    168 #endif
     165        endif !ok_gradsfile
    169166
    170167         print*
    … …  
    211208c           physic=.false.
    212209         endif
    213 
     210c correction de rappel dans couche limite
     211c F.Codron
     212         if (guide_BL) then
     213           alpha_pcor(:)=1.
     214         else
     215           do l=1,llm
     216             alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.
     217           enddo
     218         endif
    214219         itau_test=1001
    215220         step_rea=1
    … …  
    219224c    itau_test    montre si l'importation a deja ete faite au rang itau
    220225c lecture d'un fichier netcdf pour determiner le nombre de niveaux
     226         if (guide_modele) then
     227           if (ncidpl.eq.-99) ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcod)
     228         else
    221229         if (guide_u) then
    222230           if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)
    … …  
    235243         endif
    236244c
     245         endif
    237246         if (ncep) then
    238247          status=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
    … …  
    241250         endif
    242251          status=NF_INQ_DIMLEN(ncidpl,rid,nlev)
    243          print *,'nlev', nlev
    244           call ncclos(ncidpl,rcod)
     252         print *,'nlev guide', nlev
     253         call ncclos(ncidpl,rcod)
    245254c   Lecture du premier etat des reanalyses.
    246255         call read_reanalyse(1,ps
    … …  
    286295      factt=dtvr*iperiod/daysec
    287296      ztau(:)=factt/max(alpha_T(:),1.e-10)
    288 #ifdef gradsfile
     297      if (ok_gradsfile) then
    289298      call wrgrads(igrads,1,aire   ,'aire      ','aire      ' )
    290299      call wrgrads(igrads,1,dxdys  ,'dxdy      ','dxdy      ' )
    … …  
    298307      call wrgrads(igrads,llm,qrea2,'Qa        ','Qa        ' )
    299308      call wrgrads(igrads,llm,q,'Q         ','Q         ' )
    300 
    301309      call wrgrads(igrads,llm,qsat,'QSAT      ','QSAT      ' )
    302 #endif
     310      endif !(ok_gradsfile) then
    303311        endif
    304312      else
    … …  
    325333            do ij=1,ip1jmp1
    326334                a=(1.-tau)*ucovrea1(ij,l)+tau*ucovrea2(ij,l)
    327                 ucov(ij,l)=(1.-alpha_u(ij))*ucov(ij,l)+alpha_u(ij)*a
     335                ucov(ij,l)=(1.-alpha_u(ij))*ucov(ij,l)+alpha_u(ij)
     336     $                     *alpha_pcor(l)*a
    328337                if (first.and.ini_anal) ucov(ij,l)=a
    329338            enddo
    … …  
    336345            do ij=1,ip1jmp1
    337346                a=(1.-tau)*tetarea1(ij,l)+tau*tetarea2(ij,l)
    338                 teta(ij,l)=(1.-alpha_T(ij))*teta(ij,l)+alpha_T(ij)*a
     347                teta(ij,l)=(1.-alpha_T(ij))*teta(ij,l)+alpha_T(ij)
     348     $                     *alpha_pcor(l)*a
    339349                if (first.and.ini_anal) teta(ij,l)=a
    340350            enddo
    … …  
    361371c   hum relative en % -> hum specif
    362372                a=qsat(ij,l)*a*0.01
    363                 q(ij,l)=(1.-alpha_Q(ij))*q(ij,l)+alpha_Q(ij)*a
     373                q(ij,l)=(1.-alpha_Q(ij))*q(ij,l)+alpha_Q(ij)
     374     $                     *alpha_pcor(l)*a
    364375                if (first.and.ini_anal) q(ij,l)=a
    365376            enddo
    … …  
    372383            do ij=1,ip1jm
    373384                a=(1.-tau)*vcovrea1(ij,l)+tau*vcovrea2(ij,l)
    374                 vcov(ij,l)=(1.-alpha_v(ij))*vcov(ij,l)+alpha_v(ij)*a
     385                vcov(ij,l)=(1.-alpha_v(ij))*vcov(ij,l)+alpha_v(ij)
     386     $                     *alpha_pcor(l)*a
    375387                if (first.and.ini_anal) vcov(ij,l)=a
    376388            enddo
    … …  
    486498              stop
    487499            endif
     500            gamma=log(0.5)/log(gamma)
     501            if (gamma4) then
     502              gamma=min(gamma,4.)
     503            endif
    488504            print*,'gamma=',gamma
    489             gamma=log(0.5)/log(gamma)
    490505         endif
    491506      endif

  • LMDZ4/trunk/libf/dyn3d/guide.h

    r524 r1046  
    1010
    1111
    12       logical guide_u,guide_v,guide_T,guide_Q,guide_P
     12      logical guide_u,guide_v,guide_T,guide_Q,guide_P,guide_modele
     13      logical guide_BL,guide_hr,gamma4
    1314
    1415      real lat_min_guide,lat_max_guide
    … …  
    2122c     data tau_min_P,tau_max_P/0.02,10./
    2223c
    23       LOGICAL ncep,ini_anal
     24      LOGICAL ncep,ini_anal,ok_invertp
    2425      integer online
    2526
    … …  
    3738     s ncep,ini_anal,
    3839     s online,
    39      s guide_u,guide_v,guide_T,guide_Q,guide_P
     40     s guide_u,guide_v,guide_T,guide_Q,guide_P,
     41     s guide_modele,ok_invertp,
     42     s guide_hr,guide_BL,gamma4

  • LMDZ4/trunk/libf/dyn3d/logic.h

    r999 r1046  
    99      COMMON/logic/ purmats,iflag_phys,forward,leapf,apphys,            &
    1010     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
    11      &  ,read_start,ok_guide,ok_strato
     11     &  ,read_start,ok_guide,ok_strato,ok_gradsfile
    1212
    1313      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
    1414     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
    15      &  ,read_start,ok_guide,ok_strato
     15     &  ,read_start,ok_guide,ok_strato,ok_gradsfile
    1616
    1717      INTEGER iflag_phys

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

    r524 r1046  
    44c******************************************************
    55      SUBROUTINE   pres2lev(varo,varn,lmo,lmn,po,pn,
    6      %                      ni,nj)
     6     %                      ni,nj,ok_invertp)
    77c
    88c interpolation lineaire pour passer
    … …  
    1010c les variables de GCM
    1111c Francois Forget (01/1995)
    12 c
    1312c MOdif remy roca 12/97 pour passer de pres2sig
     13c Modif F.Codron 07/08 po en 3D
    1414c**********************************************************
    1515
    … …  
    2121c  ARGUMENTS
    2222c  """""""""
    23 
     23       LOGICAL ok_invertp
    2424       INTEGER lmo ! dimensions ancienne couches (input)
    2525       INTEGER lmn ! dimensions nouvelle couches (input)
    … …  
    2929       parameter(lmomx=10000,lmnmx=10000)
    3030
    31         real po(lmo)! niveau de pression en millibars
    32         real pn(ni,nj,lmn) ! niveau de pression en pascals
     31        real po(ni,nj,lmo)! niveau de pression ancienne grille
     32        real pn(ni,nj,lmn) ! niveau de pression nouvelle grille
    3333
    3434       INTEGER i,j,Nhoriz,ni,nj ! nombre de point horizontale (input)
    … …  
    4848        do i=1,ni
    4949        do j=1,nj
    50 c a chaque point de grille correspond un nouveau sigma old
    51 c qui vaut pres(l)/ps(i,j)
     50! Inversion de l'ordre des niveaux verticaux
     51          IF (ok_invertp) THEN
    5252           do lo=1,lmo
    53               zpo(lo)=po(lmo+1-lo)
     53              zpo(lo)=po(i,j,lmo+1-lo)
    5454              zvaro(lo)=varo(i,j,lmo+1-lo)
    5555           enddo
    56        
     56          ELSE
     57           do lo=1,lmo
     58              zpo(lo)=po(i,j,lo)
     59              zvaro(lo)=varo(i,j,lo)
     60           enddo
     61          ENDIF
     62
    5763           do ln=1,lmn
    5864              if (pn(i,j,ln).ge.zpo(1))then

  • 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.