Changeset 1046 for LMDZ4/trunk/libf


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
Files:
14 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

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

    r1000 r1046  
    720720      ok_strato=.FALSE.
    721721      CALL getin('ok_strato',ok_strato)
     722
     723!Config  Key  = ok_gradsfile
     724!Config  Desc = activation des sorties grads du guidage
     725!Config  Def  = n
     726!Config  Help = active les sorties grads du guidage
     727
     728       ok_gradsfile = .FALSE.
     729       CALL getin('ok_gradsfile',ok_gradsfile)
    722730
    723731      write(lunout,*)' #########################################'
    … …  
    757765      write(lunout,*)' omp_chunk = ', omp_chunk
    758766      write(lunout,*)' ok_strato = ', ok_strato
     767      write(lunout,*)' ok_gradsfile = ', ok_gradsfile
    759768c
    760769      RETURN

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

    r774 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/dyn3dpar/guide.h

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

    r985 r1046  
    7979      real alpha_T(ip1jmp1),alpha_P(ip1jmp1)
    8080      real alpha_u(ip1jmp1),alpha_v(ip1jm)
     81      real alpha_pcor(llm)
    8182      real dday_step,toto,reste,itau_test
    8283      INTEGER step_rea,count_no_rea
    … …  
    112113
    113114      save alpha_T,alpha_q,alpha_u,alpha_v,alpha_P,itau_test
     115      save alpha_pcor
    114116      save step_rea,count_no_rea
    115117
    … …  
    159161         
    160162         if (mpi_rank==0) then
     163         if (ok_gradsfile) then
    161164         call inigrads(igrads,iip1
    162165     s  ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
    163166     s  ,llm,presnivs,1.
    164167     s  ,dtgrads,file,'dyn_zon ')
     168         endif !ok_gradsfile
    165169         endif
    166170         
    … …  
    208212c           physic=.false.
    209213         endif
    210 
     214c correction de rappel dans couche limite
     215c F.Codron
     216         if (guide_BL) then
     217           alpha_pcor(:)=1.
     218         else
     219           do l=1,llm
     220             alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.
     221           enddo
     222         endif
    211223         itau_test=1001
    212224         step_rea=1
    … …  
    216228c lecture d'un fichier netcdf pour determiner le nombre de niveaux
    217229         IF (mpi_rank==0) THEN
    218        
     230         if (guide_modele) then
     231           if (ncidpl.eq.-99) ncidpl=NCOPN('apbp.nc',NCNOWRIT,rcod)
     232         else   
    219233         if (guide_u) then
    220234           if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)
    … …  
    233247         endif
    234248c
     249         endif
    235250         if (ncep) then
    236251          status=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
    … …  
    239254         endif
    240255          status=NF_INQ_DIMLEN(ncidpl,rid,nlev)
    241          print *,'nlev', nlev
    242           call ncclos(ncidpl,rcod)
    243          
    244          ENDIF
    245          
     256         print *,'nlev guide', nlev
     257         call ncclos(ncidpl,rcod)
    246258c   Lecture du premier etat des reanalyses.
    247259         call Gather_Field(ps,ip1jmp1,1,0)
    … …  
    316328     
    317329       if (mpi_rank==0) then
     330        if (ok_gradsfile) then
    318331         call wrgrads(igrads,1,aire   ,'aire      ','aire      ' )
    319332         call wrgrads(igrads,1,dxdys  ,'dxdy      ','dxdy      ' )
    … …  
    328341         call wrgrads(igrads,llm,q,'Q         ','Q         ' )
    329342         call wrgrads(igrads,llm,qsat,'QSAT      ','QSAT      ' )
     343       endif !(ok_gradsfile) then
    330344      endif
    331345     
    … …  
    357371            do ij=ijb,ije
    358372                a=(1.-tau)*ucovrea1(ij,l)+tau*ucovrea2(ij,l)
    359                 ucov(ij,l)=(1.-alpha_u(ij))*ucov(ij,l)+alpha_u(ij)*a
     373                ucov(ij,l)=(1.-alpha_u(ij))*ucov(ij,l)+alpha_u(ij)
     374     $                     *alpha_pcor(l)*a
    360375                if (first.and.ini_anal) ucov(ij,l)=a
    361376            enddo
    … …  
    368383            do ij=ijb,ije
    369384                a=(1.-tau)*tetarea1(ij,l)+tau*tetarea2(ij,l)
    370                 teta(ij,l)=(1.-alpha_T(ij))*teta(ij,l)+alpha_T(ij)*a
     385                teta(ij,l)=(1.-alpha_T(ij))*teta(ij,l)+alpha_T(ij)
     386     $                     *alpha_pcor(l)*a
    371387                if (first.and.ini_anal) teta(ij,l)=a
    372388            enddo
    … …  
    393409c   hum relative en % -> hum specif
    394410                a=qsat(ij,l)*a*0.01
    395                 q(ij,l)=(1.-alpha_Q(ij))*q(ij,l)+alpha_Q(ij)*a
     411                q(ij,l)=(1.-alpha_Q(ij))*q(ij,l)+alpha_Q(ij)
     412     $                     *alpha_pcor(l)*a
    396413                if (first.and.ini_anal) q(ij,l)=a
    397414            enddo
    … …  
    406423            do ij=ijb,ije
    407424                a=(1.-tau)*vcovrea1(ij,l)+tau*vcovrea2(ij,l)
    408                 vcov(ij,l)=(1.-alpha_v(ij))*vcov(ij,l)+alpha_v(ij)*a
     425                vcov(ij,l)=(1.-alpha_v(ij))*vcov(ij,l)+alpha_v(ij)
     426     $                     *alpha_pcor(l)*a
    409427                if (first.and.ini_anal) vcov(ij,l)=a
    410428            enddo
    … …  
    493511         enddo
    494512
    495          call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL  ')
    496          call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U     ')
     513c         call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL  ')
     514c         call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U     ')
    497515         call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v     ')
    498516
    … …  
    520538              stop
    521539            endif
     540            gamma=log(0.5)/log(gamma)
     541            if (gamma4) then
     542              gamma=min(gamma,4.)
     543            endif
    522544            print*,'gamma=',gamma
    523             gamma=log(0.5)/log(gamma)
    524545         endif
    525546      endif

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

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

    r774 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/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.