Changeset 109


Ignore:
Timestamp:
Apr 14, 2011, 11:47:04 AM (14 years ago)
Author:
slebonnois
Message:

SLebonnois: discretisation verticale: cohabitation entre
la methode Terre et les autres.

Location:
trunk
Files:
3 added
1 deleted
15 edited
5 copied
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/chantiers/commit_importants.log

    r108 r109  
    802802  Ajout dissip_horiz.tex (et .pdf) dans la documentation.
    803803
    804 Reste la discretisation verticale a voir.
    805 
     804*********************
     805**** commit_v109 ****
     806*********************
     807
     808** Discretisation verticale
     809**-----------------------
     810
     811* Modif de dyn3d[par]/iniconst.F
     812  Appel de disvert_[no]terre mis sous flag planet_type
     813
     814* Modif de dyn3d[par]/disvert.F90
     815  Renomme disvert_terre.F90
     816
     817* Ajout de dyn3d[par]/disvert_noterre.F
     818  Correspond au disvert pour Titan et Venus (et Mars ? A verifier)
     819  Declaration aps et bps enlevee (=> comvert.h)
     820
     821* Modif de dyn3d[par]/comvert.h
     822  Ajout de aps,bps
     823   => entraine la modif de limy.F(+par), vlsplt[_p].F, vlspltqs[_p].F,
     824
     825* Modif de dyn3d[par]/limy.F, vlsplt[_p].F, vlspltqs[_p].F,
     826  Changement de aps en apps (et du coup, apn en appn)
     827
     828* Ajout dyn3d[par]/exner_milieu.F (vient directement de mars/libf/dyn3d/)
     829  pour remplacer exner_hyb dans cas noterre
     830
     831* Ajout dyn3dpar/exner_milieu_p.F
     832  adapte de exner_milieu.F
     833  pour remplacer exner_hyb_p dans cas noterre
     834
     835* Modif de dyn3d[par]/leapfrog[_p].F
     836  Appels à exner_hyb[_p]/exner_milieu[_p] mis sous flag planet_type.
     837
     838* Ajout de disvert.tex[/.pdf] dans la doc
     839

  • trunk/chantiers/meschantiers-Seb.txt

    r4 r109  
    1010   Reflechir au controle des frequences de sorties
    1111
    12 - dynamique: introduire les modifs Titan/Venus dans dyn3d[par]
    13     Quid de la discretisation verticale ? Comment rendre souple ?
     12- dynamique:
     13   Il reste a voir tout ce qui est lie a l'etat initial
     14   et au newstart
    1415
    1516- I/O: evolution des I/O de Venus et Titan vers Terre ?
    1617       souplesse pour I/O Mars ?
    1718
    18 - newstart...
    19 
    2019- testphys1d...
    2120
    22 - couche eponge...
    23 

  • trunk/documentation/disvert.tex

    r108 r109  
    3333\vspace{1cm}
    3434\Large
    35 The upper boundary sponge layer
     35The vertical discretization
    3636}
    3737
    … …  
    4343\end{center}
    4444
    45 %\section{Theoretical aspects}
     45
     46\section{Theoretical aspects}
     47
     48The position of the layers:
     49\begin{itemize}
     50\item pressure limit between two layers,
     51\item pressure within the layers
     52\end{itemize}
     53
     54The Exner function: definition.
     55It corresponds to the pressure levels within the layers.
     56Used for the computation of the potential temperature.
     57For the Earth, we use a specific scheme that computes these positions so that
     58it maintains a condition of proportionality between total,
     59internal and potential energy (cf. a note from F. Hourdin).
    4660
    4761\section{Pratical aspects in the code}
    4862
    49 The sponge layer is applied at the upper boundary when the \textsf{ok\_strato}
    50 flag is set to {\em True} in \textsf{gcm.def}
    51 (this parameter also controls the application of a second step in the
    52 horizontal dissipation).
     63\begin{itemize}
     64\item \textsf{disvert\_[no]terre.F[90]}:
     65position of the interface pressure levels from an input file
     66(several possibilities).
     67Definition of ap, bp and presnivs.
     68In the planetary version, definition of aps and bps.
    5369
    54 The tendencies for the upper boundary sponge layer are computed separately in
    55 the \textsf{top\_bound.F} routine, called from \textsf{leapfrog.F}.
    56 These tendencies are \textsf{dutop}, \textsf{dvtop} and \textsf{dhtop}, in
    57 unit/s.
     70This is done only once, called at the beginning from \textsf{iniconst.F}.
    5871
    59 Three parameters may be adjusted in the \textsf{gcm.def} file:
    60 \begin{itemize}
    61 \item \textsf{iflag\_top\_bound}: selects the affected layers.
    62   \begin{itemize}
    63   \item 1: only the top 4 layers are affected. In this case, the damping rate
    64   is divided by 2 in the second layer, 4 in the third and 8 in the fourth.
    65   \item 2: layers with pressure lower than 100 times the top pressure.
    66   In this case, the damping rate depends linearly on the pressure.
    67   \end{itemize}
    68 \item \textsf{mode\_top\_bound}: selects how the fields are affected.
    69   \begin{itemize}
    70   \item 0: No sponge layer is applied.
    71   \item 1: Zonal and meridional winds are damped to zero.
    72   \item 2: Zonal and meridional winds are damped to their zonally averaged value.
    73   \item 3: Temperature, zonal and meridional winds are damped to their zonally
    74   averaged value.
    75   \end{itemize}
    76 \item \textsf{tau\_top\_bound}: damping rate (in /s) in the top layer.
     72\item Interface pressures:
     73computed in \textsf{caldyn0.F, caldyn.F, integrd.F, leapfrog.F}
     74through the \textsf{pression.F} routine.
     75
     76\item Exner function (and therefore pressure within the layers):
     77computed at three different places in \textsf{leapfrog.F} through the
     78\textsf{exner\_[hyb/milieu].F} routine.
     79For the Earth, we use \textsf{exner\_hyb.F}, that computes the positions in a
     80specific way to maintain a condition of proportionality between total,
     81internal and potential energy (cf. a note from F. Hourdin).
     82For other planets, we use \textsf{exner\_milieu.F}, that computes the positions
     83of these pressure levels exactly in the middle of each layer.
     84Though this fails to maintain the previous condition, there is no evidence of
     85any significant influence on the results, and it makes it a lot easier to
     86define correctly the level positions with the input file.
    7787\end{itemize}
    7888

  • trunk/libf/dyn3d/comvert.h

    r1 r109  
    66
    77      COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),      &
    8      &               pa,preff,nivsigs(llm),nivsig(llm+1)
     8     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
     9     &               aps(llm),bps(llm)
    910
    10       REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig
     11      REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig,aps,bps
    1112
    1213 !-----------------------------------------------------------------------

  • trunk/libf/dyn3d/disvert_terre.F90

    r107 r109  
    11! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $
    22
    3 SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
     3SUBROUTINE disvert_terre(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
    44
    55  ! Auteur : P. Le Van

  • trunk/libf/dyn3d/exner_milieu.F

    r107 r109  
    1       SUBROUTINE  exner_hyb ( ngrid, ps, p,beta, pks, pk, pkf )
     1      SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
    22c
    33c     Auteurs :  F. Forget , Y. Wanherdrick
    … …  
    1717c
    1818c     WARNING : CECI est une version speciale de exner_hyb originale
    19 c               Utilis‰ dans la version martienne pour pouvoir
    20 c               tourner avec des coordonn‰es verticales complexe
    21 c              => Il ne verifie PAS la condition la proportionalit‰ en
    22 c              ‰nergie totale/ interne / potentielle (F.Forget 2001)
     19c               Utilise dans la version martienne pour pouvoir
     20c               tourner avec des coordonnees verticales complexe
     21c              => Il ne verifie PAS la condition la proportionalite en
     22c              energie totale/ interne / potentielle (F.Forget 2001)
    2323c    ( voir note de Fr.Hourdin )  ,
    2424c

  • trunk/libf/dyn3d/iniconst.F

    r1 r109  
    5353c-----------------------------------------------------------------------
    5454
    55        CALL disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
    56 c
     55      if (planet_type.eq."earth") then
     56       CALL disvert_terre(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
     57      else
     58       CALL disvert_noterre
     59      endif
    5760c
    5861       RETURN

  • trunk/libf/dyn3d/leapfrog.F

    r108 r109  
    238238      dq(:,:,:)=0.
    239239      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    240       CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     240      if (planet_type.eq."earth") then
     241        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     242      else
     243        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     244      endif
     245
    241246c------------------
    242247c TEST PK MONOTONE
    … …  
    404409
    405410         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
    406          CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
     411         if (planet_type.eq."earth") then
     412           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     413         else
     414           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     415         endif
    407416
    408417!           rdaym_ini  = itau * dtvr / daysec
    … …  
    519528
    520529        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
    521         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     530        if (planet_type.eq."earth") then
     531          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     532        else
     533          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     534        endif
    522535
    523536

  • trunk/libf/dyn3d/limy.F

    r1 r109  
    4040      REAL qbyv(ip1jm,llm)
    4141
    42       REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2
     42      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2
    4343      Logical extremum,first
    4444      save first
    … …  
    117117
    118118c     print*,dyqv(iip1+1)
    119 c     apn=abs(dyq(1)/dyqv(iip1+1))
     119c     appn=abs(dyq(1)/dyqv(iip1+1))
    120120c     print*,dyq(ip1jm+1)
    121121c     print*,dyqv(ip1jm-iip1+1)
    122 c     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     122c     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    123123c     do ij=2,iim
    124 c        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    125 c        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
     124c        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     125c        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    126126c     enddo
    127 c     apn=min(pente_max/apn,1.)
    128 c     aps=min(pente_max/aps,1.)
     127c     appn=min(pente_max/appn,1.)
     128c     apps=min(pente_max/apps,1.)
    129129
    130130
    … …  
    132132
    133133c     if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    134 c    &   apn=0.
     134c    &   appn=0.
    135135c     if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    136136c    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    137 c    &   aps=0.
     137c    &   apps=0.
    138138
    139139c   limitation des pentes aux poles
    140140c     do ij=1,iip1
    141 c        dyq(ij)=apn*dyq(ij)
    142 c        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     141c        dyq(ij)=appn*dyq(ij)
     142c        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    143143c     enddo
    144144

  • trunk/libf/dyn3d/vlsplt.F

    r1 r109  
    478478      REAL qbyv(ip1jm,llm)
    479479
    480       REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
     480      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
    481481c     REAL newq,oldmasse
    482482      Logical extremum,first,testcpu
    … …  
    602602C     PRINT*,dyq(1)
    603603C     PRINT*,dyqv(iip1+1)
    604 C     apn=abs(dyq(1)/dyqv(iip1+1))
     604C     appn=abs(dyq(1)/dyqv(iip1+1))
    605605C     PRINT*,dyq(ip1jm+1)
    606606C     PRINT*,dyqv(ip1jm-iip1+1)
    607 C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     607C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    608608C     DO ij=2,iim
    609 C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    610 C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
     609C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     610C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    611611C     ENDDO
    612 C     apn=min(pente_max/apn,1.)
    613 C     aps=min(pente_max/aps,1.)
     612C     appn=min(pente_max/appn,1.)
     613C     apps=min(pente_max/apps,1.)
    614614C
    615615C
    … …  
    617617C
    618618C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    619 C    &   apn=0.
     619C    &   appn=0.
    620620C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    621621C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    622 C    &   aps=0.
     622C    &   apps=0.
    623623C
    624624C   limitation des pentes aux poles
    625625C     DO ij=1,iip1
    626 C        dyq(ij)=apn*dyq(ij)
    627 C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     626C        dyq(ij)=appn*dyq(ij)
     627C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    628628C     ENDDO
    629629C

  • trunk/libf/dyn3d/vlspltqs.F

    r5 r109  
    635635C     PRINT*,dyq(1)
    636636C     PRINT*,dyqv(iip1+1)
    637 C     apn=abs(dyq(1)/dyqv(iip1+1))
     637C     appn=abs(dyq(1)/dyqv(iip1+1))
    638638C     PRINT*,dyq(ip1jm+1)
    639639C     PRINT*,dyqv(ip1jm-iip1+1)
    640 C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     640C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    641641C     DO ij=2,iim
    642 C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    643 C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
     642C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     643C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    644644C     ENDDO
    645 C     apn=min(pente_max/apn,1.)
    646 C     aps=min(pente_max/aps,1.)
     645C     appn=min(pente_max/appn,1.)
     646C     apps=min(pente_max/apps,1.)
    647647C
    648648C
    … …  
    650650C
    651651C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    652 C    &   apn=0.
     652C    &   appn=0.
    653653C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    654654C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    655 C    &   aps=0.
     655C    &   apps=0.
    656656C
    657657C   limitation des pentes aux poles
    658658C     DO ij=1,iip1
    659 C        dyq(ij)=apn*dyq(ij)
    660 C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     659C        dyq(ij)=appn*dyq(ij)
     660C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    661661C     ENDDO
    662662C

  • trunk/libf/dyn3dpar/comvert.h

    r1 r109  
    11!
    2 ! $Header$
     2! $Id: comvert.h 1279 2009-12-10 09:02:56Z fairhead $
    33!
    44!-----------------------------------------------------------------------
    55!   INCLUDE 'comvert.h'
    66
    7       COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),       &
    8      &               pa,preff,nivsigs(llm),nivsig(llm+1)
     7      COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),      &
     8     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
     9     &               aps(llm),bps(llm)
    910
    10       REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig
     11      REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig,aps,bps
    1112
    12 !-----------------------------------------------------------------------
     13 !-----------------------------------------------------------------------

  • trunk/libf/dyn3dpar/disvert_terre.F90

    r107 r109  
    11! $Id: disvert.F90 1480 2011-01-31 21:29:58Z jghattas $
    22
    3 SUBROUTINE disvert(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
     3SUBROUTINE disvert_terre(pa, preff, ap, bp, dpres, presnivs, nivsigs, nivsig)
    44
    55  ! Auteur : P. Le Van

  • trunk/libf/dyn3dpar/exner_milieu.F

    r107 r109  
    1       SUBROUTINE  exner_hyb ( ngrid, ps, p,beta, pks, pk, pkf )
     1      SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
    22c
    33c     Auteurs :  F. Forget , Y. Wanherdrick
    … …  
    1717c
    1818c     WARNING : CECI est une version speciale de exner_hyb originale
    19 c               Utilis‰ dans la version martienne pour pouvoir
    20 c               tourner avec des coordonn‰es verticales complexe
    21 c              => Il ne verifie PAS la condition la proportionalit‰ en
    22 c              ‰nergie totale/ interne / potentielle (F.Forget 2001)
     19c               Utilise dans la version martienne pour pouvoir
     20c               tourner avec des coordonnees verticales complexe
     21c              => Il ne verifie PAS la condition la proportionalite en
     22c              energie totale/ interne / potentielle (F.Forget 2001)
    2323c    ( voir note de Fr.Hourdin )  ,
    2424c

  • trunk/libf/dyn3dpar/exner_milieu_p.F

    r107 r109  
    1       SUBROUTINE  exner_hyb ( ngrid, ps, p,beta, pks, pk, pkf )
     1      SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
    22c
    33c     Auteurs :  F. Forget , Y. Wanherdrick
    … …  
    1717c
    1818c     WARNING : CECI est une version speciale de exner_hyb originale
    19 c               Utilis‰ dans la version martienne pour pouvoir
    20 c               tourner avec des coordonn‰es verticales complexe
    21 c              => Il ne verifie PAS la condition la proportionalit‰ en
    22 c              ‰nergie totale/ interne / potentielle (F.Forget 2001)
     19c               Utilise dans la version martienne pour pouvoir
     20c               tourner avec des coordonnees verticales complexe
     21c              => Il ne verifie PAS la condition la proportionalite en
     22c              energie totale/ interne / potentielle (F.Forget 2001)
    2323c    ( voir note de Fr.Hourdin )  ,
    2424c
     25      USE parallel
    2526      IMPLICIT NONE
    2627c
    … …  
    4445      REAL xpn, xps
    4546      REAL SSUM
    46       EXTERNAL filtreg, SSUM
     47      EXTERNAL SSUM
     48      INTEGER ije,ijb,jje,jjb
    4749     
     50c$OMP BARRIER
     51
    4852c     -------------
    4953c     Calcul de pks
    5054c     -------------
    5155   
    52       DO   ij  = 1, ngrid
     56      ijb=ij_begin
     57      ije=ij_end
     58
     59c$OMP DO SCHEDULE(STATIC)
     60      DO   ij  = ijb, ije
    5361        pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
    5462      ENDDO
     63c$OMP ENDDO
     64c Synchro OPENMP ici
    5565
    56       DO  ij   = 1, iim
    57         ppn(ij) = aire(   ij   ) * pks(  ij     )
    58         pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    59       ENDDO
    60       xpn      = SSUM(iim,ppn,1) /apoln
    61       xps      = SSUM(iim,pps,1) /apols
    62 
    63       DO ij   = 1, iip1
    64         pks(   ij     )  =  xpn
    65         pks( ij+ip1jm )  =  xps
    66       ENDDO
     66c$OMP MASTER
     67      if (pole_nord) then
     68        DO  ij   = 1, iim
     69          ppn(ij) = aire(   ij   ) * pks(  ij     )
     70        ENDDO
     71        xpn      = SSUM(iim,ppn,1) /apoln
     72 
     73        DO ij   = 1, iip1
     74          pks(   ij     )  =  xpn
     75        ENDDO
     76      endif
     77     
     78      if (pole_sud) then
     79        DO  ij   = 1, iim
     80          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
     81        ENDDO
     82        xps      = SSUM(iim,pps,1) /apols
     83 
     84        DO ij   = 1, iip1
     85          pks( ij+ip1jm )  =  xps
     86        ENDDO
     87      endif
     88c$OMP END MASTER
    6789c
    6890c
    … …  
    7294      dum1 = cpp * (2*preff)**(-kappa)
    7395      DO l = 1, llm-1
    74         DO   ij   = 1, ngrid
     96c$OMP DO SCHEDULE(STATIC)
     97        DO   ij   = ijb, ije
    7598         pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
    7699        ENDDO
     100c$OMP ENDDO NOWAIT       
    77101      ENDDO
    78102
    … …  
    81105c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
    82106
    83       DO   ij   = 1, ngrid
     107c$OMP DO SCHEDULE(STATIC)
     108      DO   ij   = ijb, ije
    84109         pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
    85110      ENDDO
     111c$OMP ENDDO NOWAIT       
    86112
    87113
    88114c    calcul de pkf
    89115c    -------------
    90       CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    91       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
     116c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
     117      DO l = 1, llm
     118c$OMP DO SCHEDULE(STATIC)
     119         DO   ij   = ijb, ije
     120           pkf(ij,l)=pk(ij,l)
     121         ENDDO
     122c$OMP ENDDO NOWAIT             
     123      ENDDO
     124
     125c$OMP BARRIER
     126     
     127      jjb=jj_begin
     128      jje=jj_end
     129      CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
    92130     
    93131c    EST-CE UTILE ?? : calcul de beta
    94132c    --------------------------------
    95133      DO l = 2, llm
    96         DO   ij   = 1, ngrid
     134c$OMP DO SCHEDULE(STATIC)
     135        DO   ij   = ijb, ije
    97136          beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
    98137        ENDDO
     138c$OMP ENDDO NOWAIT             
    99139      ENDDO
    100140

  • trunk/libf/dyn3dpar/iniconst.F

    r1 r109  
    5353c-----------------------------------------------------------------------
    5454
    55        CALL disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
    56 c
     55      if (planet_type.eq."earth") then
     56       CALL disvert_terre(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
     57      else
     58       CALL disvert_noterre
     59      endif
    5760c
    5861       RETURN

  • trunk/libf/dyn3dpar/inidissip.F

    r1 r109  
    3131      INTEGER l,ij,idum,ii
    3232      REAL tetamin
    33       REAL pseudoz
     33      REAL Pup
    3434
    3535      REAL ran1
    … …  
    177177c   --------------------------------------------------
    178178
    179       if (ok_strato .and. llm==39) then
    180          do l=1,llm
    181             pseudoz=8.*log(preff/presnivs(l))
    182             zvert(l)=1+
    183      s      (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2.
    184      s      *(dissip_factz-1.)
    185          enddo
    186       else
    187          DO l=1,llm
    188             zvert(l)=1.
    189          ENDDO
    190          fact=2.
    191          DO l = 1, llm
    192             zz      = 1. - preff/presnivs(l)
    193             zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
    194          ENDDO
     179c First step: going from 1 to dissip_fac_mid (in gcm.def)
     180c============
     181      DO l=1,llm
     182         zz      = 1. - preff/presnivs(l)
     183         zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
     184      ENDDO
     185
     186         write(*,*) 'Dissipation : '
     187         write(*,*) 'Multiplication de la dissipation en altitude :',
     188         write(*,*) '  dissip_fac_mid =', dissip_fac_mid
     189
     190c Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
     191c==========================
     192c Utilisation de la fonction tangente hyperbolique pour augmenter
     193c arbitrairement la dissipation et donc la stabilite du modele a
     194c partir d'une certaine altitude.
     195
     196c   Le facteur multiplicatif de basse atmosphere etant deja pris
     197c   en compte, il faut diviser le facteur multiplicatif de haute
     198c   atmosphere par celui-ci.
     199
     200      if (ok_strato) then
     201
     202       Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
     203       do l=1,llm
     204         zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1)
     205     &   *(1-(0.5*(1+tanh(-6./dissip_deltaz*
     206     &    (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
     207       enddo
     208
     209         write(*,*) '  dissip_fac_up =', dissip_fac_up
     210         write(*,*) 'Transition mid /up:  Pupstart,delta =',
     211     &             dissip_pupstart,'Pa', dissip_deltaz , '(km)'
     212
    195213      endif
    196214

  • trunk/libf/dyn3dpar/leapfrog_p.F

    r108 r109  
    271271
    272272      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    273       CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     273      if (planet_type.eq."earth") then
     274        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     275      else
     276        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     277      endif
    274278c$OMP END MASTER
    275279c-----------------------------------------------------------------------
    … …  
    746750
    747751c$OMP BARRIER
    748          CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
     752         if (planet_type.eq."earth") then
     753          CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     754         else
     755          CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
     756         endif
    749757c$OMP BARRIER
    750758           jD_cur = jD_ref + day_ini - day_ref
    … …  
    10941102        CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
    10951103c$OMP BARRIER
    1096         CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     1104        if (planet_type.eq."earth") then
     1105          CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     1106        else
     1107          CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
     1108        endif
    10971109c$OMP BARRIER
    10981110

  • trunk/libf/dyn3dpar/limy.F

    r1 r109  
    4040      REAL qbyv(ip1jm,llm)
    4141
    42       REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2
     42      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2
    4343      Logical extremum,first
    4444      save first
    … …  
    116116
    117117c     print*,dyqv(iip1+1)
    118 c     apn=abs(dyq(1)/dyqv(iip1+1))
     118c     appn=abs(dyq(1)/dyqv(iip1+1))
    119119c     print*,dyq(ip1jm+1)
    120120c     print*,dyqv(ip1jm-iip1+1)
    121 c     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     121c     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    122122c     do ij=2,iim
    123 c        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    124 c        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
     123c        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     124c        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    125125c     enddo
    126 c     apn=min(pente_max/apn,1.)
    127 c     aps=min(pente_max/aps,1.)
     126c     appn=min(pente_max/appn,1.)
     127c     apps=min(pente_max/apps,1.)
    128128
    129129
    … …  
    131131
    132132c     if(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    133 c    &   apn=0.
     133c    &   appn=0.
    134134c     if(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    135135c    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    136 c    &   aps=0.
     136c    &   apps=0.
    137137
    138138c   limitation des pentes aux poles
    139139c     do ij=1,iip1
    140 c        dyq(ij)=apn*dyq(ij)
    141 c        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     140c        dyq(ij)=appn*dyq(ij)
     141c        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    142142c     enddo
    143143

  • trunk/libf/dyn3dpar/vlsplt_p.F

    r1 r109  
    561561      REAL qbyv(ip1jm,llm)
    562562
    563       REAL qpns,qpsn,apn,aps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
     563      REAL qpns,qpsn,appn,apps,dyn1,dys1,dyn2,dys2,newmasse,fn,fs
    564564c     REAL newq,oldmasse
    565565      Logical extremum,first,testcpu
    … …  
    732732C     PRINT*,dyq(1)
    733733C     PRINT*,dyqv(iip1+1)
    734 C     apn=abs(dyq(1)/dyqv(iip1+1))
     734C     appn=abs(dyq(1)/dyqv(iip1+1))
    735735C     PRINT*,dyq(ip1jm+1)
    736736C     PRINT*,dyqv(ip1jm-iip1+1)
    737 C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     737C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    738738C     DO ij=2,iim
    739 C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    740 C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
     739C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     740C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    741741C     ENDDO
    742 C     apn=min(pente_max/apn,1.)
    743 C     aps=min(pente_max/aps,1.)
     742C     appn=min(pente_max/appn,1.)
     743C     apps=min(pente_max/apps,1.)
    744744C
    745745C
    … …  
    747747C
    748748C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    749 C    &   apn=0.
     749C    &   appn=0.
    750750C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    751751C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    752 C    &   aps=0.
     752C    &   apps=0.
    753753C
    754754C   limitation des pentes aux poles
    755755C     DO ij=1,iip1
    756 C        dyq(ij)=apn*dyq(ij)
    757 C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     756C        dyq(ij)=appn*dyq(ij)
     757C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    758758C     ENDDO
    759759C

  • trunk/libf/dyn3dpar/vlspltqs_p.F

    r8 r109  
    780780C     PRINT*,dyq(1)
    781781C     PRINT*,dyqv(iip1+1)
    782 C     apn=abs(dyq(1)/dyqv(iip1+1))
     782C     appn=abs(dyq(1)/dyqv(iip1+1))
    783783C     PRINT*,dyq(ip1jm+1)
    784784C     PRINT*,dyqv(ip1jm-iip1+1)
    785 C     aps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
     785C     apps=abs(dyq(ip1jm+1)/dyqv(ip1jm-iip1+1))
    786786C     DO ij=2,iim
    787 C        apn=amax1(abs(dyq(ij)/dyqv(ij)),apn)
    788 C        aps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),aps)
     787C        appn=amax1(abs(dyq(ij)/dyqv(ij)),appn)
     788C        apps=amax1(abs(dyq(ip1jm+ij)/dyqv(ip1jm-iip1+ij)),apps)
    789789C     ENDDO
    790 C     apn=min(pente_max/apn,1.)
    791 C     aps=min(pente_max/aps,1.)
     790C     appn=min(pente_max/appn,1.)
     791C     apps=min(pente_max/apps,1.)
    792792C
    793793C
    … …  
    795795C
    796796C     IF(dyqv(ismin(iim,dyqv,1))*dyqv(ismax(iim,dyqv,1)).le.0.)
    797 C    &   apn=0.
     797C    &   appn=0.
    798798C     IF(dyqv(ismax(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1)*
    799799C    &   dyqv(ismin(iim,dyqv(ip1jm-iip1+1),1)+ip1jm-iip1+1).le.0.)
    800 C    &   aps=0.
     800C    &   apps=0.
    801801C
    802802C   limitation des pentes aux poles
    803803C     DO ij=1,iip1
    804 C        dyq(ij)=apn*dyq(ij)
    805 C        dyq(ip1jm+ij)=aps*dyq(ip1jm+ij)
     804C        dyq(ij)=appn*dyq(ij)
     805C        dyq(ip1jm+ij)=apps*dyq(ip1jm+ij)
    806806C     ENDDO
    807807C
Note: See TracChangeset for help on using the changeset viewer.