Changeset 2021


Ignore:
Timestamp:
Apr 25, 2014, 12:20:14 PM (10 years ago)
Author:
lguez
Message:

Removed unused variables pks, pk, pkf from main program unit gcm.

Encapsulated procedures exner_hyb, exner_hyb_p, exner_hyb_loc,
exner_milieu, exner_milieu_p and exner_milieu_loc into
modules. (Compulsory to allow optional arguments.)

In the procedures exner_hyb, exner_hyb_p, exner_hyb_loc, donwgraded
arguments alpha and beta to local variables. In exner_milieu,
exner_milieu_p and exner_milieu_loc, removed beta altogether. In the
six procedures exner_*, made pkf an optional argument. Made some
cosmetic modifications in order to keep the six procedures exner_* as
close as possible.

In the six procedures exner_*, removed the averaging of pks at the
poles: this is not useful because ps is already the same at all
longitudes at the poles. This modification changes the results of the
program. Motivation: had to do this for exner_hyb because we call it
from test_disvert with a few surface pressure values.

In all the procedures calling exner_*, removed the variables alpha and
beta. Also removed variables alpha and beta from module leapfrog_mod
and from module call_calfis_mod.

Removed actual argument pkf in call to exner_hyb* and exner_milieu*
from guide_interp, guide_main, iniacademic and iniacademic_loc (pkf
was not used in those procedures).

Argument workvar of startget_dyn is used only if varname is tpot or

  1. When varname is tpot or q, the actual argument associated to

workvar in etat0_netcdf is not y. So y in etat0_netcdf is a
place-holder, never used. So we remove optional argument y in the
calls to exner_hyb and exner_milieu from etat0_netcdf.

Created procedure test_disvert, called only by etat0_netcdf. This
procedure tests the order of pressure values at half-levels and full
levels.

Location:
LMDZ5/trunk/libf
Files:
1 added
15 edited
6 moved

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3d/gcm.F

    r1930 r2021  
    105105      REAL ps(ip1jmp1)                       ! pression  au sol
    106106      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
    107       REAL pks(ip1jmp1)                      ! exner au  sol
    108       REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    109       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    110107      REAL masse(ip1jmp1,llm)                ! masse d'air
    111108      REAL phis(ip1jmp1)                     ! geopotentiel au sol
    … …  
    131128      data call_iniphys/.true./
    132129
    133       REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
    134130c+jld variables test conservation energie
    135131c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)

  • LMDZ5/trunk/libf/dyn3d/guide_mod.F90

    r1907 r2021  
    589589  SUBROUTINE guide_interp(psi,teta)
    590590 
     591  use exner_hyb_m, only: exner_hyb
     592  use exner_milieu_m, only: exner_milieu
    591593  IMPLICIT NONE
    592594
    … …  
    610612  REAL, DIMENSION (iip1,jjm,llm)     :: pbary
    611613  ! Variables pour fonction Exner (P milieu couche)
    612   REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
    613   REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
     614  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
    614615  REAL, DIMENSION (iip1,jjp1)        :: pks   
    615616  REAL                               :: prefkap,unskap
    … …  
    676677    CALL pression( ip1jmp1, ap, bp, psi, p )
    677678    if (pressure_exner) then
    678       CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
     679      CALL exner_hyb(ip1jmp1,psi,p,pks,pk)
    679680    else
    680       CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf)
     681      CALL exner_milieu(ip1jmp1,psi,p,pks,pk)
    681682    endif
    682683!    ....  Calcul de pls , pression au milieu des couches ,en Pascals

  • LMDZ5/trunk/libf/dyn3d/iniacademic.F90

    r1907 r2021  
    1414#endif
    1515  USE Write_Field
     16  use exner_hyb_m, only: exner_hyb
     17  use exner_milieu_m, only: exner_milieu
    1618
    1719  !   Author:    Frederic Hourdin      original: 15/01/93
    … …  
    5456  REAL pks(ip1jmp1)                      ! exner au  sol
    5557  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    56   REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    5758  REAL phi(ip1jmp1,llm)                  ! geopotentiel
    5859  REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
    … …  
    7071  integer idum
    7172
    72   REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
     73  REAL zdtvr
    7374 
    7475  character(len=*),parameter :: modname="iniacademic"
    … …  
    223224        CALL pression ( ip1jmp1, ap, bp, ps, p       )
    224225        if (pressure_exner) then
    225           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    226         else
    227           call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
     226          CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
     227        else
     228          call exner_milieu(ip1jmp1,ps,p,pks,pk)
    228229        endif
    229230        CALL massdair(p,masse)

  • LMDZ5/trunk/libf/dyn3d/leapfrog.F

    r1987 r2021  
    1919     &                       iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins,
    2020     &                       periodav, ok_dyn_ave, output_grads_dyn
     21      use exner_hyb_m, only: exner_hyb
     22      use exner_milieu_m, only: exner_milieu
     23
    2124      IMPLICIT NONE
    2225
    … …  
    158161      character*10 string10
    159162
    160       REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
    161163      REAL :: flxw(ip1jmp1,llm)  ! flux de masse verticale
    162164
    … …  
    217219      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    218220      if (pressure_exner) then
    219         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     221        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
    220222      else
    221         CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     223        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
    222224      endif
    223225
    … …  
    373375         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
    374376         if (pressure_exner) then
    375            CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
     377           CALL exner_hyb(  ip1jmp1, ps, p,pks, pk, pkf )
    376378         else
    377            CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     379           CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
    378380         endif
    379381
    … …  
    448450          CALL massdair(p,masse)
    449451          if (pressure_exner) then
    450             CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
     452            CALL exner_hyb(ip1jmp1,ps,p,pks,pk,pkf)
    451453          else
    452             CALL exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
     454            CALL exner_milieu(ip1jmp1,ps,p,pks,pk,pkf)
    453455          endif
    454456
    … …  
    506508        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
    507509        if (pressure_exner) then
    508           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     510          CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
    509511        else
    510           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     512          CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
    511513        endif
    512514        CALL massdair(p,masse)

  • LMDZ5/trunk/libf/dyn3d_common/exner_hyb_m.F90

    r1992 r2021  
    1 !
    2 ! $Id $
    3 !
    4       SUBROUTINE  exner_hyb ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
    5 c
    6 c     Auteurs :  P.Le Van  , Fr. Hourdin  .
    7 c    ..........
    8 c
    9 c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
    10 c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
    11 c
    12 c   ************************************************************************
    13 c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des
    14 c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
    15 c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
    16 c   ************************************************************************
    17 c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
    18 c    la pression et la fonction d'Exner  au  sol  .
    19 c
    20 c                                 -------- z                                   
    21 c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
    22 c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
    23 c    ( voir note de Fr.Hourdin )  ,
    24 c
    25 c    on determine successivement , du haut vers le bas des couches, les
    26 c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
    27 c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
    28 c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
    29 c
    30 c
    31       IMPLICIT NONE
    32 c
    33 #include "dimensions.h"
    34 #include "paramet.h"
    35 #include "comconst.h"
    36 #include "comgeom.h"
    37 #include "comvert.h"
    38 #include "serre.h"
     1module exner_hyb_m
    392
    40       INTEGER  ngrid
    41       REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
    42       REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
     3  IMPLICIT NONE
    434
    44 c    .... variables locales   ...
     5contains
    456
    46       INTEGER l, ij
    47       REAL unpl2k,dellta
     7  SUBROUTINE  exner_hyb ( ngrid, ps, p, pks, pk, pkf )
    488
    49       REAL ppn(iim),pps(iim)
    50       REAL xpn, xps
    51       REAL SSUM
    52 c
    53       logical,save :: firstcall=.true.
    54       character(len=*),parameter :: modname="exner_hyb"
    55      
    56       ! Sanity check
    57       if (firstcall) then
    58         ! sanity checks for Shallow Water case (1 vertical layer)
    59         if (llm.eq.1) then
     9    !     Auteurs :  P.Le Van  , Fr. Hourdin  .
     10    !    ..........
     11    !
     12    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
     13    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
     14    !
     15    !   ************************************************************************
     16    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
     17    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
     18    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
     19    !   ************************************************************************
     20    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
     21    !    la pression et la fonction d'Exner  au  sol  .
     22    !
     23    !                                 -------- z
     24    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
     25    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
     26    !    ( voir note de Fr.Hourdin )  ,
     27    !
     28    !    on determine successivement , du haut vers le bas des couches, les
     29    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
     30    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
     31    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
     32    !
     33    !
     34    !
     35    include "dimensions.h"
     36    include "paramet.h"
     37    include "comconst.h"
     38    include "comgeom.h"
     39    include "comvert.h"
     40    include "serre.h"
     41
     42    INTEGER  ngrid
     43    REAL p(ngrid,llmp1),pk(ngrid,llm)
     44    real, optional:: pkf(ngrid,llm)
     45    REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
     46
     47    !    .... variables locales   ...
     48
     49    INTEGER l, ij
     50    REAL unpl2k,dellta
     51
     52    logical,save :: firstcall=.true.
     53    character(len=*),parameter :: modname="exner_hyb"
     54
     55    ! Sanity check
     56    if (firstcall) then
     57       ! sanity checks for Shallow Water case (1 vertical layer)
     58       if (llm.eq.1) then
    6059          if (kappa.ne.1) then
    61             call abort_gcm(modname,
    62      &      "kappa!=1 , but running in Shallow Water mode!!",42)
     60             call abort_gcm(modname, &
     61                  "kappa!=1 , but running in Shallow Water mode!!",42)
    6362          endif
    6463          if (cpp.ne.r) then
    65             call abort_gcm(modname,
    66      &      "cpp!=r , but running in Shallow Water mode!!",42)
     64             call abort_gcm(modname, &
     65                  "cpp!=r , but running in Shallow Water mode!!",42)
    6766          endif
    68         endif ! of if (llm.eq.1)
     67       endif ! of if (llm.eq.1)
    6968
    70         firstcall=.false.
    71       endif ! of if (firstcall)
     69       firstcall=.false.
     70    endif ! of if (firstcall)
    7271
    73       if (llm.eq.1) then
    74        
    75         ! Compute pks(:),pk(:),pkf(:)
    76        
    77         DO   ij  = 1, ngrid
    78           pks(ij) = (cpp/preff) * ps(ij)
     72    ! Specific behaviour for Shallow Water (1 vertical layer) case:
     73    if (llm.eq.1) then
     74
     75       ! Compute pks(:),pk(:),pkf(:)
     76
     77       DO   ij  = 1, ngrid
     78          pks(ij) = (cpp/preff) * ps(ij)
    7979          pk(ij,1) = .5*pks(ij)
    80         ENDDO
    81        
    82         CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    83         CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
    84        
    85         ! our work is done, exit routine
    86         return
     80       ENDDO
    8781
    88       endif ! of if (llm.eq.1)
     82       if (present(pkf)) then
     83          pkf = pk
     84          CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
     85       end if
    8986
    90 !!!! General case:
    91      
    92       unpl2k    = 1.+ 2.* kappa
    93 c
    94       DO   ij  = 1, ngrid
    95         pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
    96       ENDDO
     87       ! our work is done, exit routine
     88       return
     89    endif ! of if (llm.eq.1)
    9790
    98       DO  ij   = 1, iim
    99         ppn(ij) = aire(   ij   ) * pks(  ij     )
    100         pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    101       ENDDO
    102       xpn      = SSUM(iim,ppn,1) /apoln
    103       xps      = SSUM(iim,pps,1) /apols
     91    ! General case:
    10492
    105       DO ij   = 1, iip1
    106         pks(   ij     )  =  xpn
    107         pks( ij+ip1jm )  =  xps
    108       ENDDO
    109 c
    110 c
    111 c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
    112 c
    113       DO     ij      = 1, ngrid
     93    unpl2k    = 1.+ 2.* kappa
     94
     95    !     -------------
     96    !     Calcul de pks
     97    !     -------------
     98
     99    DO   ij  = 1, ngrid
     100       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
     101    ENDDO
     102
     103    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
     104    !
     105    DO     ij      = 1, ngrid
    114106       alpha(ij,llm) = 0.
    115107       beta (ij,llm) = 1./ unpl2k
    116       ENDDO
    117 c
    118 c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
    119 c
    120       DO l = llm -1 , 2 , -1
    121 c
    122         DO ij = 1, ngrid
    123         dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
    124         alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
    125         beta (ij,l)  =   p(ij,l  ) / dellta   
    126         ENDDO
    127 c
    128       ENDDO
    129 c
    130 c  ***********************************************************************
    131 c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
    132 c
    133       DO   ij   = 1, ngrid
    134        pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
    135      *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
    136       ENDDO
    137 c
    138 c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
    139 c
    140       DO l = 2, llm
    141         DO   ij   = 1, ngrid
    142          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
    143         ENDDO
    144       ENDDO
    145 c
    146 c
    147       CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    148       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
    149      
     108    ENDDO
     109    !
     110    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
     111    !
     112    DO l = llm -1 , 2 , -1
     113       !
     114       DO ij = 1, ngrid
     115          dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
     116          alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
     117          beta (ij,l)  =   p(ij,l  ) / dellta   
     118       ENDDO
     119    ENDDO
    150120
    151       RETURN
    152       END
     121    !  ***********************************************************************
     122    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
     123    !
     124    DO   ij   = 1, ngrid
     125       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  / &
     126            (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
     127    ENDDO
     128    !
     129    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
     130    !
     131    DO l = 2, llm
     132       DO   ij   = 1, ngrid
     133          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
     134       ENDDO
     135    ENDDO
     136
     137    if (present(pkf)) then
     138       !    calcul de pkf
     139       pkf = pk
     140       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
     141    end if
     142
     143  END SUBROUTINE exner_hyb
     144
     145end module exner_hyb_m

  • LMDZ5/trunk/libf/dyn3d_common/exner_milieu_m.F90

    r1992 r2021  
    1 !
    2 ! $Id $
    3 !
    4       SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
    5 c
    6 c     Auteurs :  F. Forget , Y. Wanherdrick
    7 c P.Le Van  , Fr. Hourdin  .
    8 c    ..........
    9 c
    10 c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
    11 c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
    12 c
    13 c   ************************************************************************
    14 c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
    15 c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
    16 c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
    17 c   ************************************************************************
    18 c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
    19 c    la pression et la fonction d'Exner  au  sol  .
    20 c
    21 c     WARNING : CECI est une version speciale de exner_hyb originale
    22 c               Utilise dans la version martienne pour pouvoir
    23 c               tourner avec des coordonnees verticales complexe
    24 c              => Il ne verifie PAS la condition la proportionalite en
    25 c              energie totale/ interne / potentielle (F.Forget 2001)
    26 c    ( voir note de Fr.Hourdin )  ,
    27 c
    28       IMPLICIT NONE
    29 c
    30 #include "dimensions.h"
    31 #include "paramet.h"
    32 #include "comconst.h"
    33 #include "comgeom.h"
    34 #include "comvert.h"
    35 #include "serre.h"
     1module exner_milieu_m
    362
    37       INTEGER  ngrid
    38       REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
    39       REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
     3  IMPLICIT NONE
    404
    41 c    .... variables locales   ...
     5contains
    426
    43       INTEGER l, ij
    44       REAL dum1
     7  SUBROUTINE  exner_milieu ( ngrid, ps, p, pks, pk, pkf )
     8    !
     9    !     Auteurs :  F. Forget , Y. Wanherdrick
     10    ! P.Le Van  , Fr. Hourdin  .
     11    !    ..........
     12    !
     13    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
     14    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
     15    !
     16    !   ************************************************************************
     17    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
     18    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
     19    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
     20    !   ************************************************************************
     21    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
     22    !    la pression et la fonction d'Exner  au  sol  .
     23    !
     24    !     WARNING : CECI est une version speciale de exner_hyb originale
     25    !               Utilise dans la version martienne pour pouvoir
     26    !               tourner avec des coordonnees verticales complexe
     27    !              => Il ne verifie PAS la condition la proportionalite en
     28    !              energie totale/ interne / potentielle (F.Forget 2001)
     29    !    ( voir note de Fr.Hourdin )  ,
     30    !
     31    !
     32    include "dimensions.h"
     33    include "paramet.h"
     34    include "comconst.h"
     35    include "comgeom.h"
     36    include "comvert.h"
     37    include "serre.h"
    4538
    46       REAL ppn(iim),pps(iim)
    47       REAL xpn, xps
    48       REAL SSUM
    49       EXTERNAL SSUM
    50       logical,save :: firstcall=.true.
    51       character(len=*),parameter :: modname="exner_milieu"
     39    INTEGER  ngrid
     40    REAL p(ngrid,llmp1),pk(ngrid,llm)
     41    real, optional:: pkf(ngrid,llm)
     42    REAL ps(ngrid),pks(ngrid)
    5243
    53       ! Sanity check
    54       if (firstcall) then
    55         ! sanity checks for Shallow Water case (1 vertical layer)
    56         if (llm.eq.1) then
     44    !    .... variables locales   ...
     45
     46    INTEGER l, ij
     47    REAL dum1
     48
     49    logical,save :: firstcall=.true.
     50    character(len=*),parameter :: modname="exner_milieu"
     51
     52    ! Sanity check
     53    if (firstcall) then
     54       ! sanity checks for Shallow Water case (1 vertical layer)
     55       if (llm.eq.1) then
    5756          if (kappa.ne.1) then
    58             call abort_gcm(modname,
    59      &      "kappa!=1 , but running in Shallow Water mode!!",42)
     57             call abort_gcm(modname, &
     58                  "kappa!=1 , but running in Shallow Water mode!!",42)
    6059          endif
    6160          if (cpp.ne.r) then
    62             call abort_gcm(modname,
    63      &      "cpp!=r , but running in Shallow Water mode!!",42)
     61             call abort_gcm(modname, &
     62                  "cpp!=r , but running in Shallow Water mode!!",42)
    6463          endif
    65         endif ! of if (llm.eq.1)
     64       endif ! of if (llm.eq.1)
    6665
    67         firstcall=.false.
    68       endif ! of if (firstcall)
     66       firstcall=.false.
     67    endif ! of if (firstcall)
    6968
    70 !!!! Specific behaviour for Shallow Water (1 vertical layer) case:
    71       if (llm.eq.1) then
    72      
    73         ! Compute pks(:),pk(:),pkf(:)
    74        
    75         DO   ij  = 1, ngrid
    76           pks(ij) = (cpp/preff) * ps(ij) 
     69    ! Specific behaviour for Shallow Water (1 vertical layer) case:
     70    if (llm.eq.1) then
     71
     72       ! Compute pks(:),pk(:),pkf(:)
     73
     74       DO   ij  = 1, ngrid
     75          pks(ij) = (cpp/preff) * ps(ij)
    7776          pk(ij,1) = .5*pks(ij)
    78         ENDDO
    79        
    80         CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    81         CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
    82        
    83         ! our work is done, exit routine
    84         return
     77       ENDDO
    8578
    86       endif ! of if (llm.eq.1)
     79       if (present(pkf)) then
     80          pkf = pk
     81          CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
     82       end if
    8783
    88 !!!! General case:
     84       ! our work is done, exit routine
     85       return
     86    endif ! of if (llm.eq.1)
    8987
    90 c     -------------
    91 c     Calcul de pks
    92 c     -------------
    93    
    94       DO   ij  = 1, ngrid
    95         pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
    96       ENDDO
     88    ! General case:
    9789
    98       DO  ij   = 1, iim
    99         ppn(ij) = aire(   ij   ) * pks(  ij     )
    100         pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    101       ENDDO
    102       xpn      = SSUM(iim,ppn,1) /apoln
    103       xps      = SSUM(iim,pps,1) /apols
     90    !     -------------
     91    !     Calcul de pks
     92    !     -------------
    10493
    105       DO ij   = 1, iip1
    106         pks(   ij     )  =  xpn
    107         pks( ij+ip1jm )  =  xps
    108       ENDDO
    109 c
    110 c
    111 c    .... Calcul de pk  pour la couche l
    112 c    --------------------------------------------
    113 c
    114       dum1 = cpp * (2*preff)**(-kappa)
    115       DO l = 1, llm-1
    116         DO   ij   = 1, ngrid
    117          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
    118         ENDDO
    119       ENDDO
     94    DO   ij  = 1, ngrid
     95       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
     96    ENDDO
    12097
    121 c    .... Calcul de pk  pour la couche l = llm ..
    122 c    (on met la meme distance (en log pression)  entre Pk(llm)
    123 c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
     98    !    .... Calcul de pk  pour la couche l
     99    !    --------------------------------------------
     100    !
     101    dum1 = cpp * (2*preff)**(-kappa)
     102    DO l = 1, llm-1
     103       DO   ij   = 1, ngrid
     104          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
     105       ENDDO
     106    ENDDO
    124107
    125       DO   ij   = 1, ngrid
    126          pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
    127       ENDDO
     108    !    .... Calcul de pk  pour la couche l = llm ..
     109    !    (on met la meme distance (en log pression)  entre Pk(llm)
     110    !    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
    128111
     112    DO   ij   = 1, ngrid
     113       pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
     114    ENDDO
    129115
    130 c    calcul de pkf
    131 c    -------------
    132       CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    133       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
    134      
    135 c    EST-CE UTILE ?? : calcul de beta
    136 c    --------------------------------
    137       DO l = 2, llm
    138         DO   ij   = 1, ngrid
    139           beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
    140         ENDDO
    141       ENDDO
     116    if (present(pkf)) then
     117       !    calcul de pkf
     118       pkf = pk
     119       CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
     120    end if
    142121
    143       RETURN
    144       END
     122  END SUBROUTINE exner_milieu
     123
     124end module exner_milieu_m

  • LMDZ5/trunk/libf/dyn3dmem/call_calfis_mod.F90

    r1987 r2021  
    1212
    1313    REAL,POINTER,SAVE :: p(:,:)
    14     REAL,POINTER,SAVE :: alpha(:,:)
    15     REAL,POINTER,SAVE :: beta(:,:)
    1614    REAL,POINTER,SAVE :: pks(:)
    1715    REAL,POINTER,SAVE :: pk(:,:)
    … …  
    5351    CALL allocate_u(flxw,llm,d)
    5452    CALL allocate_u(p,llmp1,d)
    55     CALL allocate_u(alpha,llm,d)
    56     CALL allocate_u(beta,llm,d)
    5753    CALL allocate_u(pks,d)
    5854    CALL allocate_u(pk,llm,d)
    … …  
    7571                         phis_dyn,q_dyn,flxw_dyn)
    7672  USE dimensions_mod
     73  use exner_hyb_loc_m, only: exner_hyb_loc
     74  use exner_milieu_loc_m, only: exner_milieu_loc
    7775  USE parallel_lmdz
    7876  USE times
    … …  
    201199
    202200  !$OMP BARRIER
    203     CALL exner_hyb_loc(  ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     201    CALL exner_hyb_loc(  ip1jmp1, ps, p, pks, pk, pkf )
    204202  !$OMP BARRIER
    205203    CALL geopot_loc  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
    … …  
    343341!$OMP BARRIER
    344342    if (pressure_exner) then
    345       CALL exner_hyb_loc(ijnb_u,ps,p,alpha,beta,pks,pk,pkf)
     343      CALL exner_hyb_loc(ijnb_u,ps,p,pks,pk,pkf)
    346344    else
    347       CALL exner_milieu_loc(ijnb_u,ps,p,beta,pks,pk,pkf)
     345      CALL exner_milieu_loc(ijnb_u,ps,p,pks,pk,pkf)
    348346    endif
    349347!$OMP BARRIER

  • LMDZ5/trunk/libf/dyn3dmem/exner_hyb_loc_m.F90

    r1992 r2021  
    1 c
    2 c $Id$
    3 c
    4       SUBROUTINE  exner_hyb_loc(ngrid, ps, p,alpha,beta, pks,pk,pkf)
    5 c
    6 c     Auteurs :  P.Le Van  , Fr. Hourdin  .
    7 c    ..........
    8 c
    9 c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
    10 c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
    11 c
    12 c   ************************************************************************
    13 c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des
    14 c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
    15 c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
    16 c   ************************************************************************
    17 c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
    18 c    la pression et la fonction d'Exner  au  sol  .
    19 c
    20 c                                 -------- z                                   
    21 c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
    22 c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
    23 c    ( voir note de Fr.Hourdin )  ,
    24 c
    25 c    on determine successivement , du haut vers le bas des couches, les
    26 c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
    27 c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
    28 c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
    29 c
    30 c
    31       USE parallel_lmdz
    32       USE mod_filtreg_p
    33       USE write_field_loc
    34       IMPLICIT NONE
    35 c
    36 #include "dimensions.h"
    37 #include "paramet.h"
    38 #include "comconst.h"
    39 #include "comgeom.h"
    40 #include "comvert.h"
    41 #include "serre.h"
     1module exner_hyb_loc_m
    422
    43       INTEGER  ngrid
    44       REAL p(ijb_u:ije_u,llmp1),pk(ijb_u:ije_u,llm)
    45       REAL pkf(ijb_u:ije_u,llm)
    46       REAL ps(ijb_u:ije_u),pks(ijb_u:ije_u)
    47       REAL alpha(ijb_u:ije_u,llm),beta(ijb_u:ije_u,llm)
     3  IMPLICIT NONE
    484
    49 c    .... variables locales   ...
     5contains
    506
    51       INTEGER l, ij
    52       REAL unpl2k,dellta
     7  SUBROUTINE  exner_hyb_loc(ngrid, ps, p, pks,pk,pkf)
    538
    54       REAL ppn(iim),pps(iim)
    55       REAL xpn, xps
    56       REAL SSUM
    57       EXTERNAL SSUM
    58       INTEGER ije,ijb,jje,jjb
    59       logical,save :: firstcall=.true.
    60 !$OMP THREADPRIVATE(firstcall)
    61       character(len=*),parameter :: modname="exner_hyb_loc"
    62 c
    63 c$OMP BARRIER           
     9    !     Auteurs :  P.Le Van  , Fr. Hourdin  .
     10    !    ..........
     11    !
     12    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
     13    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
     14    !
     15    !   ************************************************************************
     16    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
     17    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
     18    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
     19    !   ************************************************************************
     20    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
     21    !    la pression et la fonction d'Exner  au  sol  .
     22    !
     23    !                                 -------- z
     24    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
     25    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
     26    !    ( voir note de Fr.Hourdin )  ,
     27    !
     28    !    on determine successivement , du haut vers le bas des couches, les
     29    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
     30    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
     31    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
     32    !
     33    !
     34    USE parallel_lmdz
     35    USE mod_filtreg_p
     36    USE write_field_loc
     37    !
     38    include "dimensions.h"
     39    include "paramet.h"
     40    include "comconst.h"
     41    include "comgeom.h"
     42    include "comvert.h"
     43    include "serre.h"
    6444
    65       ! Sanity check
    66       if (firstcall) then
    67         ! sanity checks for Shallow Water case (1 vertical layer)
    68         if (llm.eq.1) then
     45    INTEGER  ngrid
     46    REAL p(ijb_u:ije_u,llmp1),pk(ijb_u:ije_u,llm)
     47    REAL, optional:: pkf(ijb_u:ije_u,llm)
     48    REAL ps(ijb_u:ije_u),pks(ijb_u:ije_u)
     49    REAL alpha(ijb_u:ije_u,llm),beta(ijb_u:ije_u,llm)
     50
     51    !    .... variables locales   ...
     52
     53    INTEGER l, ij
     54    REAL unpl2k,dellta
     55
     56    INTEGER ije,ijb,jje,jjb
     57    logical,save :: firstcall=.true.
     58    !$OMP THREADPRIVATE(firstcall)
     59    character(len=*),parameter :: modname="exner_hyb_loc"
     60    !
     61    !$OMP BARRIER           
     62
     63    ! Sanity check
     64    if (firstcall) then
     65       ! sanity checks for Shallow Water case (1 vertical layer)
     66       if (llm.eq.1) then
    6967          if (kappa.ne.1) then
    70             call abort_gcm(modname,
    71      &      "kappa!=1 , but running in Shallow Water mode!!",42)
     68             call abort_gcm(modname, &
     69                  "kappa!=1 , but running in Shallow Water mode!!",42)
    7270          endif
    7371          if (cpp.ne.r) then
    74             call abort_gcm(modname,
    75      &      "cpp!=r , but running in Shallow Water mode!!",42)
     72             call abort_gcm(modname, &
     73                  "cpp!=r , but running in Shallow Water mode!!",42)
    7674          endif
    77         endif ! of if (llm.eq.1)
     75       endif ! of if (llm.eq.1)
    7876
    79         firstcall=.false.
    80       endif ! of if (firstcall)
     77       firstcall=.false.
     78    endif ! of if (firstcall)
    8179
    82 c$OMP BARRIER
     80    !$OMP BARRIER
    8381
    84 ! Specific behaviour for Shallow Water (1 vertical layer) case
    85       if (llm.eq.1) then
    86        
    87         ! Compute pks(:),pk(:),pkf(:)
    88         ijb=ij_begin
    89         ije=ij_end
    90 !$OMP DO SCHEDULE(STATIC)
    91         DO ij=ijb, ije
    92           pks(ij)=(cpp/preff)*ps(ij)
     82    ! Specific behaviour for Shallow Water (1 vertical layer) case:
     83    if (llm.eq.1) then
     84
     85       ! Compute pks(:),pk(:),pkf(:)
     86       ijb=ij_begin
     87       ije=ij_end
     88       !$OMP DO SCHEDULE(STATIC)
     89       DO ij=ijb, ije
     90          pks(ij) = (cpp/preff) * ps(ij)
    9391          pk(ij,1) = .5*pks(ij)
    94           pkf(ij,1)=pk(ij,1)
    95         ENDDO
    96 !$OMP ENDDO
     92          if (present(pkf)) pkf(ij,1)=pk(ij,1)
     93       ENDDO
     94       !$OMP ENDDO
    9795
    98 !$OMP MASTER
    99       if (pole_nord) then
    100         DO  ij   = 1, iim
    101           ppn(ij) = aire(   ij   ) * pks(  ij     )
    102         ENDDO
    103         xpn      = SSUM(iim,ppn,1) /apoln
    104  
    105         DO ij   = 1, iip1
    106           pks(   ij     )  =  xpn
    107           pk(ij,1) = .5*pks(ij)
    108           pkf(ij,1)=pk(ij,1)
    109         ENDDO
    110       endif
    111      
    112       if (pole_sud) then
    113         DO  ij   = 1, iim
    114           pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    115         ENDDO
    116         xps      = SSUM(iim,pps,1) /apols
    117  
    118         DO ij   = 1, iip1
    119           pks( ij+ip1jm )  =  xps
    120           pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
    121           pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
    122         ENDDO
    123       endif
    124 !$OMP END MASTER
    125 !$OMP BARRIER
    126         jjb=jj_begin
    127         jje=jj_end
    128         CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm,
    129      &                 2, 1, .TRUE., 1 )
     96       !$OMP BARRIER
     97       if (present(pkf)) then
     98          jjb=jj_begin
     99          jje=jj_end
     100          CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, &
     101               2, 1, .TRUE., 1 )
     102       end if
    130103
    131         ! our work is done, exit routine
    132         return
    133       endif ! of if (llm.eq.1)
     104       ! our work is done, exit routine
     105       return
     106    endif ! of if (llm.eq.1)
    134107
     108    ! General case:
    135109
    136       unpl2k    = 1.+ 2.* kappa
    137 c
    138       ijb=ij_begin
    139       ije=ij_end
     110    unpl2k    = 1.+ 2.* kappa
    140111
    141 c$OMP DO SCHEDULE(STATIC)
    142       DO   ij  = ijb, ije
    143         pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
    144       ENDDO
    145 c$OMP ENDDO
    146 c Synchro OPENMP ici
     112    !     -------------
     113    !     Calcul de pks
     114    !     -------------
    147115
    148 c$OMP MASTER
    149       if (pole_nord) then
    150         DO  ij   = 1, iim
    151           ppn(ij) = aire(   ij   ) * pks(  ij     )
    152         ENDDO
    153         xpn      = SSUM(iim,ppn,1) /apoln
    154  
    155         DO ij   = 1, iip1
    156           pks(   ij     )  =  xpn
    157         ENDDO
    158       endif
    159      
    160       if (pole_sud) then
    161         DO  ij   = 1, iim
    162           pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    163         ENDDO
    164         xps      = SSUM(iim,pps,1) /apols
    165  
    166         DO ij   = 1, iip1
    167           pks( ij+ip1jm )  =  xps
    168         ENDDO
    169       endif
    170 c$OMP END MASTER
    171 c$OMP BARRIER
    172 c
    173 c
    174 c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
    175 c
    176 c$OMP DO SCHEDULE(STATIC)
    177       DO     ij      = ijb,ije
     116    ijb=ij_begin
     117    ije=ij_end
     118
     119    !$OMP DO SCHEDULE(STATIC)
     120    DO   ij  = ijb, ije
     121       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
     122    ENDDO
     123    !$OMP ENDDO
     124    ! Synchro OPENMP ici
     125
     126    !$OMP BARRIER
     127    !
     128    !
     129    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
     130    !
     131    !$OMP DO SCHEDULE(STATIC)
     132    DO     ij      = ijb,ije
    178133       alpha(ij,llm) = 0.
    179134       beta (ij,llm) = 1./ unpl2k
    180       ENDDO
    181 c$OMP ENDDO NOWAIT
    182 c
    183 c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
    184 c
    185       DO l = llm -1 , 2 , -1
    186 c
    187 c$OMP DO SCHEDULE(STATIC)
    188         DO ij = ijb, ije
    189         dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
    190         alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
    191         beta (ij,l)  =   p(ij,l  ) / dellta   
    192         ENDDO
    193 c$OMP ENDDO NOWAIT
    194 c
    195       ENDDO
     135    ENDDO
     136    !$OMP ENDDO NOWAIT
     137    !
     138    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
     139    !
     140    DO l = llm -1 , 2 , -1
     141       !
     142       !$OMP DO SCHEDULE(STATIC)
     143       DO ij = ijb, ije
     144          dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
     145          alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
     146          beta (ij,l)  =   p(ij,l  ) / dellta   
     147       ENDDO
     148       !$OMP ENDDO NOWAIT
     149    ENDDO
    196150
    197 c
    198 c  ***********************************************************************
    199 c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
    200 c
    201 c$OMP DO SCHEDULE(STATIC)
    202       DO   ij   = ijb, ije
    203        pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
    204      *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
    205       ENDDO
    206 c$OMP ENDDO NOWAIT
    207 c
    208 c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
    209 c
    210       DO l = 2, llm
    211 c$OMP DO SCHEDULE(STATIC)
    212         DO   ij   = ijb, ije
    213          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
    214         ENDDO
    215 c$OMP ENDDO NOWAIT       
    216       ENDDO
    217 c
    218 c
    219 c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    220       DO l = 1, llm
    221 c$OMP DO SCHEDULE(STATIC)
    222          DO   ij   = ijb, ije
    223            pkf(ij,l)=pk(ij,l)
    224          ENDDO
    225 c$OMP ENDDO NOWAIT             
    226       ENDDO
     151    !  ***********************************************************************
     152    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
     153    !
     154    !$OMP DO SCHEDULE(STATIC)
     155    DO   ij   = ijb, ije
     156       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  / &
     157            (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
     158    ENDDO
     159    !$OMP ENDDO NOWAIT
     160    !
     161    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
     162    !
     163    DO l = 2, llm
     164       !$OMP DO SCHEDULE(STATIC)
     165       DO   ij   = ijb, ije
     166          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
     167       ENDDO
     168       !$OMP ENDDO NOWAIT       
     169    ENDDO
    227170
    228 c$OMP BARRIER
    229      
    230       jjb=jj_begin
    231       jje=jj_end
     171    if (present(pkf)) then
     172       !    calcul de pkf
     173
     174       DO l = 1, llm
     175          !$OMP DO SCHEDULE(STATIC)
     176          DO   ij   = ijb, ije
     177             pkf(ij,l)=pk(ij,l)
     178          ENDDO
     179          !$OMP ENDDO NOWAIT             
     180       ENDDO
     181
     182       !$OMP BARRIER
     183
     184       jjb=jj_begin
     185       jje=jj_end
    232186#ifdef DEBUG_IO   
    233       call WriteField_u('pkf',pkf)
     187       call WriteField_u('pkf',pkf)
    234188#endif
    235       CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm,
    236      &                 2, 1, .TRUE., 1 )
     189       CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, &
     190            2, 1, .TRUE., 1 )
    237191#ifdef DEBUG_IO   
    238       call WriteField_u('pkf',pkf)
    239 #endif     
     192       call WriteField_u('pkf',pkf)
     193#endif     
     194    end if
    240195
    241       RETURN
    242       END
     196  END SUBROUTINE exner_hyb_loc
     197
     198end module exner_hyb_loc_m

  • LMDZ5/trunk/libf/dyn3dmem/exner_milieu_loc_m.F90

    r1992 r2021  
    1 !
    2 ! $Id$
    3 !
    4       SUBROUTINE  exner_milieu_loc ( ngrid, ps, p,beta, pks, pk, pkf )
    5 c
    6 c     Auteurs :  F. Forget , Y. Wanherdrick
    7 c P.Le Van  , Fr. Hourdin  .
    8 c    ..........
    9 c
    10 c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
    11 c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
    12 c
    13 c   ************************************************************************
    14 c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
    15 c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
    16 c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
    17 c   ************************************************************************
    18 c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
    19 c    la pression et la fonction d'Exner  au  sol  .
    20 c
    21 c     WARNING : CECI est une version speciale de exner_hyb originale
    22 c               Utilise dans la version martienne pour pouvoir
    23 c               tourner avec des coordonnees verticales complexe
    24 c              => Il ne verifie PAS la condition la proportionalite en
    25 c              energie totale/ interne / potentielle (F.Forget 2001)
    26 c    ( voir note de Fr.Hourdin )  ,
    27 c
    28       USE parallel_lmdz
    29       USE mod_filtreg_p
    30       IMPLICIT NONE
    31 c
    32 #include "dimensions.h"
    33 #include "paramet.h"
    34 #include "comconst.h"
    35 #include "comgeom.h"
    36 #include "comvert.h"
    37 #include "serre.h"
     1module exner_milieu_loc_m
    382
    39       INTEGER  ngrid
    40       REAL p(ijb_u:ije_u,llmp1),pk(ijb_u:ije_u,llm)
    41       REAL pkf(ijb_u:ije_u,llm)
    42       REAL ps(ijb_u:ije_u),pks(ijb_u:ije_u)
    43       REAL alpha(ijb_u:ije_u,llm),beta(ijb_u:ije_u,llm)
     3  IMPLICIT NONE
    444
    45 c    .... variables locales   ...
     5contains
    466
    47       INTEGER l, ij
    48       REAL dum1
     7  SUBROUTINE  exner_milieu_loc ( ngrid, ps, p, pks, pk, pkf )
     8    !
     9    !     Auteurs :  F. Forget , Y. Wanherdrick
     10    ! P.Le Van  , Fr. Hourdin  .
     11    !    ..........
     12    !
     13    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
     14    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
     15    !
     16    !   ************************************************************************
     17    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
     18    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
     19    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
     20    !   ************************************************************************
     21    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
     22    !    la pression et la fonction d'Exner  au  sol  .
     23    !
     24    !     WARNING : CECI est une version speciale de exner_hyb originale
     25    !               Utilise dans la version martienne pour pouvoir
     26    !               tourner avec des coordonnees verticales complexe
     27    !              => Il ne verifie PAS la condition la proportionalite en
     28    !              energie totale/ interne / potentielle (F.Forget 2001)
     29    !    ( voir note de Fr.Hourdin )  ,
     30    !
     31    USE parallel_lmdz
     32    USE mod_filtreg_p
     33    !
     34    include "dimensions.h"
     35    include "paramet.h"
     36    include "comconst.h"
     37    include "comgeom.h"
     38    include "comvert.h"
     39    include "serre.h"
    4940
    50       REAL ppn(iim),pps(iim)
    51       REAL xpn, xps
    52       REAL SSUM
    53       EXTERNAL SSUM
    54       INTEGER ije,ijb,jje,jjb
    55       logical,save :: firstcall=.true.
    56 !$OMP THREADPRIVATE(firstcall)
    57       character(len=*),parameter :: modname="exner_milieu_loc"
     41    INTEGER  ngrid
     42    REAL p(ijb_u:ije_u,llmp1),pk(ijb_u:ije_u,llm)
     43    REAL, optional:: pkf(ijb_u:ije_u,llm)
     44    REAL ps(ijb_u:ije_u),pks(ijb_u:ije_u)
    5845
    59       ! Sanity check
    60       if (firstcall) then
    61        
    62         ! sanity checks for Shallow Water case (1 vertical layer)
    63         if (llm.eq.1) then
     46    !    .... variables locales   ...
     47
     48    INTEGER l, ij
     49    REAL dum1
     50
     51    INTEGER ije,ijb,jje,jjb
     52    logical,save :: firstcall=.true.
     53    !$OMP THREADPRIVATE(firstcall)
     54    character(len=*),parameter :: modname="exner_milieu_loc"
     55
     56    ! Sanity check
     57    if (firstcall) then
     58       ! sanity checks for Shallow Water case (1 vertical layer)
     59       if (llm.eq.1) then
    6460          if (kappa.ne.1) then
    65             call abort_gcm(modname,
    66      &      "kappa!=1 , but running in Shallow Water mode!!",42)
     61             call abort_gcm(modname, &
     62                  "kappa!=1 , but running in Shallow Water mode!!",42)
    6763          endif
    6864          if (cpp.ne.r) then
    69             call abort_gcm(modname,
    70      &      "cpp!=r , but running in Shallow Water mode!!",42)
     65             call abort_gcm(modname, &
     66                  "cpp!=r , but running in Shallow Water mode!!",42)
    7167          endif
    72         endif ! of if (llm.eq.1)
     68       endif ! of if (llm.eq.1)
    7369
    74         firstcall=.false.
    75       endif ! of if (firstcall)
    76      
    77 c$OMP BARRIER
     70       firstcall=.false.
     71    endif ! of if (firstcall)
    7872
    79 ! Specific behaviour for Shallow Water (1 vertical layer) case
    80       if (llm.eq.1) then
    81              
    82         ! Compute pks(:),pk(:),pkf(:)
    83         ijb=ij_begin
    84         ije=ij_end
    85 !$OMP DO SCHEDULE(STATIC)
    86         DO ij=ijb, ije
    87           pks(ij)=(cpp/preff)*ps(ij)
     73    !$OMP BARRIER
     74
     75    ! Specific behaviour for Shallow Water (1 vertical layer) case:
     76    if (llm.eq.1) then
     77
     78       ! Compute pks(:),pk(:),pkf(:)
     79       ijb=ij_begin
     80       ije=ij_end
     81       !$OMP DO SCHEDULE(STATIC)
     82       DO ij=ijb, ije
     83          pks(ij) = (cpp/preff) * ps(ij)
    8884          pk(ij,1) = .5*pks(ij)
    89           pkf(ij,1)=pk(ij,1)
    90         ENDDO
    91 !$OMP ENDDO
     85          if (present(pkf)) pkf(ij,1)=pk(ij,1)
     86       ENDDO
     87       !$OMP ENDDO
    9288
    93 !$OMP MASTER
    94       if (pole_nord) then
    95         DO  ij   = 1, iim
    96           ppn(ij) = aire(   ij   ) * pks(  ij     )
    97         ENDDO
    98         xpn      = SSUM(iim,ppn,1) /apoln
    99  
    100         DO ij   = 1, iip1
    101           pks(   ij     )  =  xpn
    102           pk(ij,1) = .5*pks(ij)
    103           pkf(ij,1)=pk(ij,1)
    104         ENDDO
    105       endif
    106      
    107       if (pole_sud) then
    108         DO  ij   = 1, iim
    109           pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    110         ENDDO
    111         xps      = SSUM(iim,pps,1) /apols
    112  
    113         DO ij   = 1, iip1
    114           pks( ij+ip1jm )  =  xps
    115           pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
    116           pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
    117         ENDDO
    118       endif
    119 !$OMP END MASTER
    120 !$OMP BARRIER
    121         jjb=jj_begin
    122         jje=jj_end
    123         CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm,
    124      &                 2, 1, .TRUE., 1 )
     89       !$OMP BARRIER
     90       if (present(pkf)) then
     91          jjb=jj_begin
     92          jje=jj_end
     93          CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, &
     94               2, 1, .TRUE., 1 )
     95       end if
    12596
    126         ! our work is done, exit routine
    127         return
    128       endif ! of if (llm.eq.1)
     97       ! our work is done, exit routine
     98       return
     99    endif ! of if (llm.eq.1)
    129100
    130 !!!! General case:
     101    ! General case:
    131102
    132 c     -------------
    133 c     Calcul de pks
    134 c     -------------
    135    
    136       ijb=ij_begin
    137       ije=ij_end
     103    !     -------------
     104    !     Calcul de pks
     105    !     -------------
    138106
    139 c$OMP DO SCHEDULE(STATIC)
    140       DO   ij  = ijb, ije
    141         pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
    142       ENDDO
    143 c$OMP ENDDO
    144 c Synchro OPENMP ici
     107    ijb=ij_begin
     108    ije=ij_end
    145109
    146 c$OMP MASTER
    147       if (pole_nord) then
    148         DO  ij   = 1, iim
    149           ppn(ij) = aire(   ij   ) * pks(  ij     )
    150         ENDDO
    151         xpn      = SSUM(iim,ppn,1) /apoln
    152  
    153         DO ij   = 1, iip1
    154           pks(   ij     )  =  xpn
    155         ENDDO
    156       endif
    157      
    158       if (pole_sud) then
    159         DO  ij   = 1, iim
    160           pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    161         ENDDO
    162         xps      = SSUM(iim,pps,1) /apols
    163  
    164         DO ij   = 1, iip1
    165           pks( ij+ip1jm )  =  xps
    166         ENDDO
    167       endif
    168 c$OMP END MASTER
    169 c$OMP BARRIER
    170 c
    171 c
    172 c    .... Calcul de pk  pour la couche l
    173 c    --------------------------------------------
    174 c
    175       dum1 = cpp * (2*preff)**(-kappa)
    176       DO l = 1, llm-1
    177 c$OMP DO SCHEDULE(STATIC)
    178         DO   ij   = ijb, ije
    179          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
    180         ENDDO
    181 c$OMP ENDDO NOWAIT       
    182       ENDDO
     110    !$OMP DO SCHEDULE(STATIC)
     111    DO   ij  = ijb, ije
     112       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
     113    ENDDO
     114    !$OMP ENDDO
     115    ! Synchro OPENMP ici
    183116
    184 c    .... Calcul de pk  pour la couche l = llm ..
    185 c    (on met la meme distance (en log pression)  entre Pk(llm)
    186 c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
     117    !$OMP BARRIER
     118    !
     119    !
     120    !    .... Calcul de pk  pour la couche l
     121    !    --------------------------------------------
     122    !
     123    dum1 = cpp * (2*preff)**(-kappa)
     124    DO l = 1, llm-1
     125       !$OMP DO SCHEDULE(STATIC)
     126       DO   ij   = ijb, ije
     127          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
     128       ENDDO
     129       !$OMP ENDDO NOWAIT
     130    ENDDO
    187131
    188 c$OMP DO SCHEDULE(STATIC)
    189       DO   ij   = ijb, ije
    190          pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
    191       ENDDO
    192 c$OMP ENDDO NOWAIT       
     132    !    .... Calcul de pk  pour la couche l = llm ..
     133    !    (on met la meme distance (en log pression)  entre Pk(llm)
     134    !    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
    193135
     136    !$OMP DO SCHEDULE(STATIC)
     137    DO   ij   = ijb, ije
     138       pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
     139    ENDDO
     140    !$OMP ENDDO NOWAIT       
    194141
    195 c    calcul de pkf
    196 c    -------------
    197 c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    198       DO l = 1, llm
    199 c$OMP DO SCHEDULE(STATIC)
    200          DO   ij   = ijb, ije
    201            pkf(ij,l)=pk(ij,l)
    202          ENDDO
    203 c$OMP ENDDO NOWAIT             
    204       ENDDO
     142    if (present(pkf)) then
     143       !    calcul de pkf
    205144
    206 c$OMP BARRIER
    207      
    208       jjb=jj_begin
    209       jje=jj_end
    210       CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm,
    211      &                 2, 1, .TRUE., 1 )
    212      
    213 c    EST-CE UTILE ?? : calcul de beta
    214 c    --------------------------------
    215       DO l = 2, llm
    216 c$OMP DO SCHEDULE(STATIC)
    217         DO   ij   = ijb, ije
    218           beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
    219         ENDDO
    220 c$OMP ENDDO NOWAIT             
    221       ENDDO
     145       DO l = 1, llm
     146          !$OMP DO SCHEDULE(STATIC)
     147          DO   ij   = ijb, ije
     148             pkf(ij,l)=pk(ij,l)
     149          ENDDO
     150          !$OMP ENDDO NOWAIT
     151       ENDDO
    222152
    223       RETURN
    224       END
     153       !$OMP BARRIER
     154
     155       jjb=jj_begin
     156       jje=jj_end
     157       CALL filtreg_p ( pkf,jjb_u,jje_u,jjb,jje, jmp1, llm, &
     158            2, 1, .TRUE., 1 )
     159    end if
     160
     161  END SUBROUTINE exner_milieu_loc
     162
     163end module exner_milieu_loc_m

  • LMDZ5/trunk/libf/dyn3dmem/gcm.F

    r1995 r2021  
    9898      REAL,ALLOCATABLE,SAVE  :: ps(:)         ! pression  au sol
    9999c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
    100 c      REAL pks(ip1jmp1)                      ! exner au  sol
    101 c      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    102 c      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    103100      REAL,ALLOCATABLE,SAVE  :: masse(:,:)    ! masse d'air
    104101      REAL,ALLOCATABLE,SAVE  :: phis(:)       ! geopotentiel au sol
    … …  
    124121      data call_iniphys/.true./
    125122
    126 c      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
    127123c+jld variables test conservation energie
    128124c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)

  • LMDZ5/trunk/libf/dyn3dmem/guide_loc_mod.F90

    r1907 r2021  
    329329!=======================================================================
    330330  SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
     331    use exner_hyb_loc_m, only: exner_hyb_loc
     332    use exner_milieu_loc_m, only: exner_milieu_loc
    331333    USE parallel_lmdz
    332334    USE control_mod
    … …  
    353355    REAL, ALLOCATABLE, SAVE, DIMENSION (:,:) :: f_addv ! var aux: champ de guidage
    354356    ! Variables pour fonction Exner (P milieu couche)
    355     REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:)    :: pk, pkf
    356     REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:)    :: alpha, beta
     357    REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:)    :: pk
    357358    REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)        :: pks   
    358359    REAL                               :: unskap
    … …  
    399400        ALLOCATE(f_addv(ijb_v:ije_v,llm) )
    400401        ALLOCATE(pk(iip1,jjb_u:jje_u,llm)  )
    401         ALLOCATE(pkf(iip1,jjb_u:jje_u,llm)  )
    402         ALLOCATE(alpha(iip1,jjb_u:jje_u,llm)  )
    403         ALLOCATE(beta(iip1,jjb_u:jje_u,llm)  )
    404402        ALLOCATE(pks(iip1,jjb_u:jje_u)  )
    405403        ALLOCATE(p(ijb_u:ije_u,llmp1) )
    … …  
    539537        CALL pression_loc( ijnb_u, ap, bp, ps, p )
    540538        if (pressure_exner) then
    541           CALL exner_hyb_loc(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
     539          CALL exner_hyb_loc(ip1jmp1,ps,p,pks,pk)
    542540        else
    543           CALL exner_milieu_loc(ip1jmp1,ps,p,beta,pks,pk,pkf)
     541          CALL exner_milieu_loc(ip1jmp1,ps,p,pks,pk)
    544542        endif
    545543!$OMP BARRIER       
    … …  
    894892!=======================================================================
    895893  SUBROUTINE guide_interp(psi,teta)
     894    use exner_hyb_loc_m, only: exner_hyb_loc
     895    use exner_milieu_loc_m, only: exner_milieu_loc
    896896  USE parallel_lmdz
    897897  USE mod_hallo
    … …  
    919919  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)     :: pbary
    920920  ! Variables pour fonction Exner (P milieu couche)
    921   REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)    :: pk, pkf
    922   REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:)    :: alpha, beta
     921  REAL, ALLOCATABLE, SAVE,DIMENSION (:,:,:)    :: pk
    923922  REAL ,ALLOCATABLE, SAVE, DIMENSION (:,:)        :: pks   
    924923  REAL                               :: unskap
    … …  
    949948      ALLOCATE(pbary(iip1,jjb_v:jje_v,llm) )   
    950949      ALLOCATE(pk(iip1,jjb_u:jje_u,llm) )   
    951       ALLOCATE(pkf(iip1,jjb_u:jje_u,llm)  )   
    952       ALLOCATE(alpha(iip1,jjb_u:jje_u,llm) )   
    953       ALLOCATE(beta(iip1,jjb_u:jje_u,llm) )   
    954950      ALLOCATE(pks (iip1,jjb_u:jje_u) )   
    955951      ALLOCATE(qsat(ijb_u:ije_u,llm) )   
    … …  
    10321028        CALL pression_loc( ijnb_u, ap, bp, psi, p )
    10331029        if (disvert_type==1) then
    1034           CALL exner_hyb_loc(ijnb_u,psi,p,alpha,beta,pks,pk,pkf)
     1030          CALL exner_hyb_loc(ijnb_u,psi,p,pks,pk)
    10351031        else ! we assume that we are in the disvert_type==2 case
    1036           CALL exner_milieu_loc(ijnb_u,psi,p,beta,pks,pk,pkf)
     1032          CALL exner_milieu_loc(ijnb_u,psi,p,pks,pk)
    10371033        endif
    10381034        unskap=1./kappa

  • LMDZ5/trunk/libf/dyn3dmem/iniacademic_loc.F90

    r1907 r2021  
    44SUBROUTINE iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0)
    55
     6  use exner_hyb_m, only: exner_hyb
     7  use exner_milieu_m, only: exner_milieu
    68  USE filtreg_mod
    79  USE infotrac, ONLY : nqtot
    … …  
    5860  REAL pks(ip1jmp1)                      ! exner au  sol
    5961  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    60   REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    6162  REAL phi(ip1jmp1,llm)                  ! geopotentiel
    6263  REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
    … …  
    7576
    7677  REAL zdtvr
    77   real,allocatable :: alpha(:,:),beta(:,:)
    7878 
    7979  character(len=*),parameter :: modname="iniacademic"
    … …  
    219219       allocate(masse_glo(ip1jmp1,llm))
    220220       allocate(phis_glo(ip1jmp1))
    221        allocate(alpha(ip1jmp1,llm))
    222        allocate(beta(ip1jmp1,llm))
    223221
    224222        ! surface pressure
    … …  
    238236        CALL pression ( ip1jmp1, ap, bp, ps_glo, p       )
    239237        if (pressure_exner) then
    240           CALL exner_hyb( ip1jmp1, ps_glo, p,alpha,beta, pks, pk, pkf )
     238          CALL exner_hyb( ip1jmp1, ps_glo, p, pks, pk )
    241239        else
    242           call exner_milieu(ip1jmp1,ps_glo,p,beta,pks,pk,pkf)
     240          call exner_milieu(ip1jmp1,ps_glo,p,pks,pk)
    243241        endif
    244242        CALL massdair(p,masse_glo)
    … …  
    301299        deallocate(ps_glo)
    302300        deallocate(phis_glo)
    303         deallocate(alpha)
    304         deallocate(beta)
    305301     ENDIF ! of IF (.NOT. read_start)
    306302  endif academic_case

  • LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F

    r1987 r2021  
    3131       USE call_calfis_mod, ONLY : call_calfis
    3232       USE leapfrog_mod
     33       use exner_hyb_loc_m, only: exner_hyb_loc
     34       use exner_milieu_loc_m, only: exner_milieu_loc
    3335      IMPLICIT NONE
    3436
    … …  
    156158      character*10 string10
    157159
    158 !      REAL,SAVE,ALLOCATABLE :: alpha(:,:),beta(:,:)
    159160!      REAL,SAVE,ALLOCATABLE :: flxw(:,:) ! flux de masse verticale
    160161
    … …  
    261262!      ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
    262263!      ALLOCATE(finvmaold(ijb_u:ije_u,llm))
    263 !      ALLOCATE(alpha(ijb_u:ije_u,llm),beta(ijb_u:ije_u,llm))
    264264!      ALLOCATE(flxw(ijb_u:ije_u,llm))
    265265!      ALLOCATE(ecin(ijb_u:ije_u,llm),ecin0(ijb_u:ije_u,llm))
    … …  
    284284c$OMP END MASTER
    285285      if (pressure_exner) then
    286       CALL exner_hyb_loc( ijnb_u, ps, p,alpha,beta, pks, pk, pkf)
     286      CALL exner_hyb_loc( ijnb_u, ps, p, pks, pk, pkf)
    287287      else
    288         CALL exner_milieu_loc( ijnb_u, ps, p, beta, pks, pk, pkf )
     288        CALL exner_milieu_loc( ijnb_u, ps, p, pks, pk, pkf )
    289289      endif
    290290c-----------------------------------------------------------------------
    … …  
    780780
    781781! c$OMP BARRIER
    782 !          CALL exner_hyb_loc(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
     782!          CALL exner_hyb_loc(  ip1jmp1, ps, p,pks, pk, pkf )
    783783! c$OMP BARRIER
    784784!            jD_cur = jD_ref + day_ini - day_ref
    … …  
    11351135c$OMP BARRIER
    11361136        if (pressure_exner) then
    1137         CALL exner_hyb_loc( ijnb_u, ps, p,alpha,beta, pks, pk, pkf )
     1137        CALL exner_hyb_loc( ijnb_u, ps, p, pks, pk, pkf )
    11381138        else
    1139           CALL exner_milieu_loc( ijnb_u, ps, p, beta, pks, pk, pkf )
     1139          CALL exner_milieu_loc( ijnb_u, ps, p, pks, pk, pkf )
    11401140        endif
    11411141c$OMP BARRIER

  • LMDZ5/trunk/libf/dyn3dmem/leapfrog_mod.F90

    r1987 r2021  
    2727  REAL,POINTER,SAVE :: dq(:,:,:)
    2828  REAL,POINTER,SAVE :: finvmaold(:,:)
    29   REAL,POINTER,SAVE :: alpha(:,:)
    30   REAL,POINTER,SAVE :: beta(:,:)
    3129  REAL,POINTER,SAVE :: flxw(:,:)
    3230  REAL,POINTER,SAVE :: unat(:,:)
    … …  
    7977    CALL allocate_u(dq,llm,nqtot,d)
    8078    CALL allocate_u(finvmaold,llm,d)
    81     CALL allocate_u(alpha,llm,d)
    82     CALL allocate_u(beta,llm,d)
    8379    CALL allocate_u(flxw,llm,d)
    8480    CALL allocate_u(unat,llm,d)
    … …  
    129125    CALL switch_u(dq,distrib_caldyn,dist)
    130126    CALL switch_u(finvmaold,distrib_caldyn,dist)
    131     CALL switch_u(alpha,distrib_caldyn,dist)
    132     CALL switch_u(beta,distrib_caldyn,dist)
    133127    CALL switch_u(flxw,distrib_caldyn,dist)
    134128    CALL switch_u(unat,distrib_caldyn,dist)

  • LMDZ5/trunk/libf/dyn3dpar/exner_hyb_p_m.F90

    r1992 r2021  
    1 !
    2 ! $Id $
    3 !
    4       SUBROUTINE  exner_hyb_p ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
    5 c
    6 c     Auteurs :  P.Le Van  , Fr. Hourdin  .
    7 c    ..........
    8 c
    9 c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
    10 c    .... alpha,beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
    11 c
    12 c   ************************************************************************
    13 c    Calcule la fonction d'Exner pk = Cp * p ** kappa , aux milieux des
    14 c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
    15 c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
    16 c   ************************************************************************
    17 c  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
    18 c    la pression et la fonction d'Exner  au  sol  .
    19 c
    20 c                                 -------- z                                   
    21 c    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
    22 c                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
    23 c    ( voir note de Fr.Hourdin )  ,
    24 c
    25 c    on determine successivement , du haut vers le bas des couches, les
    26 c    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
    27 c    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
    28 c     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
    29 c
    30 c
    31       USE parallel_lmdz
    32       IMPLICIT NONE
    33 c
    34 #include "dimensions.h"
    35 #include "paramet.h"
    36 #include "comconst.h"
    37 #include "comgeom.h"
    38 #include "comvert.h"
    39 #include "serre.h"
     1module exner_hyb_p_m
    402
    41       INTEGER  ngrid
    42       REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
    43       REAL ps(ngrid),pks(ngrid), alpha(ngrid,llm),beta(ngrid,llm)
     3  IMPLICIT NONE
    444
    45 c    .... variables locales   ...
     5contains
    466
    47       INTEGER l, ij
    48       REAL unpl2k,dellta
     7  SUBROUTINE  exner_hyb_p ( ngrid, ps, p, pks, pk, pkf )
    498
    50       REAL ppn(iim),pps(iim)
    51       REAL xpn, xps
    52       REAL SSUM
    53       EXTERNAL SSUM
    54       INTEGER ije,ijb,jje,jjb
    55       logical,save :: firstcall=.true.
    56 !$OMP THREADPRIVATE(firstcall)
    57       character(len=*),parameter :: modname="exner_hyb_p"
    58 c
     9    !     Auteurs :  P.Le Van  , Fr. Hourdin  .
     10    !    ..........
     11    !
     12    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
     13    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
     14    !
     15    !   ************************************************************************
     16    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
     17    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
     18    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
     19    !   ************************************************************************
     20    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
     21    !    la pression et la fonction d'Exner  au  sol  .
     22    !
     23    !                                 -------- z
     24    !    A partir des relations  ( 1 ) p*dz(pk) = kappa *pk*dz(p)      et
     25    !                            ( 2 ) pk(l) = alpha(l)+ beta(l)*pk(l-1)
     26    !    ( voir note de Fr.Hourdin )  ,
     27    !
     28    !    on determine successivement , du haut vers le bas des couches, les
     29    !    coef. alpha(llm),beta(llm) .,.,alpha(l),beta(l),,,alpha(2),beta(2),
     30    !    puis pk(ij,1). Ensuite ,on calcule,du bas vers le haut des couches, 
     31    !     pk(ij,l)  donne  par la relation (2),  pour l = 2 a l = llm .
     32    !
     33    !
     34    USE parallel_lmdz
     35    !
     36    include "dimensions.h"
     37    include "paramet.h"
     38    include "comconst.h"
     39    include "comgeom.h"
     40    include "comvert.h"
     41    include "serre.h"
    5942
    60       ! Sanity check
    61       if (firstcall) then
    62         ! sanity checks for Shallow Water case (1 vertical layer)
    63         if (llm.eq.1) then
     43    INTEGER  ngrid
     44    REAL p(ngrid,llmp1),pk(ngrid,llm)
     45    REAL, optional:: pkf(ngrid,llm)
     46    REAL ps(ngrid),pks(ngrid)
     47    REAL alpha(ngrid,llm),beta(ngrid,llm)
     48
     49    !    .... variables locales   ...
     50
     51    INTEGER l, ij
     52    REAL unpl2k,dellta
     53
     54    INTEGER ije,ijb,jje,jjb
     55    logical,save :: firstcall=.true.
     56    !$OMP THREADPRIVATE(firstcall)
     57    character(len=*),parameter :: modname="exner_hyb_p"
     58
     59    ! Sanity check
     60    if (firstcall) then
     61       ! sanity checks for Shallow Water case (1 vertical layer)
     62       if (llm.eq.1) then
    6463          if (kappa.ne.1) then
    65             call abort_gcm(modname,
    66      &      "kappa!=1 , but running in Shallow Water mode!!",42)
     64             call abort_gcm(modname, &
     65                  "kappa!=1 , but running in Shallow Water mode!!",42)
    6766          endif
    6867          if (cpp.ne.r) then
    69             call abort_gcm(modname,
    70      &      "cpp!=r , but running in Shallow Water mode!!",42)
     68             call abort_gcm(modname, &
     69                  "cpp!=r , but running in Shallow Water mode!!",42)
    7170          endif
    72         endif ! of if (llm.eq.1)
     71       endif ! of if (llm.eq.1)
    7372
    74         firstcall=.false.
    75       endif ! of if (firstcall)
     73       firstcall=.false.
     74    endif ! of if (firstcall)
    7675
    77 c$OMP BARRIER
     76    !$OMP BARRIER
    7877
    79 ! Specific behaviour for Shallow Water (1 vertical layer) case
    80       if (llm.eq.1) then
    81      
    82         ! Compute pks(:),pk(:),pkf(:)
    83         ijb=ij_begin
    84         ije=ij_end
    85 !$OMP DO SCHEDULE(STATIC)
    86         DO ij=ijb, ije
    87           pks(ij)=(cpp/preff)*ps(ij)
     78    ! Specific behaviour for Shallow Water (1 vertical layer) case:
     79    if (llm.eq.1) then
     80
     81       ! Compute pks(:),pk(:),pkf(:)
     82       ijb=ij_begin
     83       ije=ij_end
     84       !$OMP DO SCHEDULE(STATIC)
     85       DO ij=ijb, ije
     86          pks(ij) = (cpp/preff) * ps(ij)
    8887          pk(ij,1) = .5*pks(ij)
    89           pkf(ij,1)=pk(ij,1)
    90         ENDDO
    91 !$OMP ENDDO
     88          if (present(pkf)) pkf(ij,1)=pk(ij,1)
     89       ENDDO
     90       !$OMP ENDDO
    9291
    93 !$OMP MASTER
    94       if (pole_nord) then
    95         DO  ij   = 1, iim
    96           ppn(ij) = aire(   ij   ) * pks(  ij     )
    97         ENDDO
    98         xpn      = SSUM(iim,ppn,1) /apoln
    99  
    100         DO ij   = 1, iip1
    101           pks(   ij     )  =  xpn
    102           pk(ij,1) = .5*pks(ij)
    103           pkf(ij,1)=pk(ij,1)
    104         ENDDO
    105       endif
    106      
    107       if (pole_sud) then
    108         DO  ij   = 1, iim
    109           pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    110         ENDDO
    111         xps      = SSUM(iim,pps,1) /apols
    112  
    113         DO ij   = 1, iip1
    114           pks( ij+ip1jm )  =  xps
    115           pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
    116           pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
    117         ENDDO
    118       endif
    119 !$OMP END MASTER
    120 !$OMP BARRIER
    121         jjb=jj_begin
    122         jje=jj_end
    123         CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
     92       !$OMP BARRIER
     93       if (present(pkf)) then
     94          jjb=jj_begin
     95          jje=jj_end
     96          CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
     97       end if
    12498
    125         ! our work is done, exit routine
    126         return
    127       endif ! of if (llm.eq.1)
     99       ! our work is done, exit routine
     100       return
     101    endif ! of if (llm.eq.1)
    128102
    129 !!!! General case:
     103    ! General case:
    130104
    131       unpl2k    = 1.+ 2.* kappa
    132 c
    133       ijb=ij_begin
    134       ije=ij_end
     105    unpl2k    = 1.+ 2.* kappa
    135106
    136 c$OMP DO SCHEDULE(STATIC)
    137       DO   ij  = ijb, ije
    138         pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
    139       ENDDO
    140 c$OMP ENDDO
    141 c Synchro OPENMP ici
     107    !     -------------
     108    !     Calcul de pks
     109    !     -------------
    142110
    143 c$OMP MASTER
    144       if (pole_nord) then
    145         DO  ij   = 1, iim
    146           ppn(ij) = aire(   ij   ) * pks(  ij     )
    147         ENDDO
    148         xpn      = SSUM(iim,ppn,1) /apoln
    149  
    150         DO ij   = 1, iip1
    151           pks(   ij     )  =  xpn
    152         ENDDO
    153       endif
    154      
    155       if (pole_sud) then
    156         DO  ij   = 1, iim
    157           pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    158         ENDDO
    159         xps      = SSUM(iim,pps,1) /apols
    160  
    161         DO ij   = 1, iip1
    162           pks( ij+ip1jm )  =  xps
    163         ENDDO
    164       endif
    165 c$OMP END MASTER
    166 c$OMP BARRIER
    167 c
    168 c
    169 c    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
    170 c
    171 c$OMP DO SCHEDULE(STATIC)
    172       DO     ij      = ijb,ije
     111    ijb=ij_begin
     112    ije=ij_end
     113
     114    !$OMP DO SCHEDULE(STATIC)
     115    DO   ij  = ijb, ije
     116       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
     117    ENDDO
     118    !$OMP ENDDO
     119    ! Synchro OPENMP ici
     120
     121    !$OMP BARRIER
     122    !
     123    !
     124    !    .... Calcul des coeff. alpha et beta  pour la couche l = llm ..
     125    !
     126    !$OMP DO SCHEDULE(STATIC)
     127    DO     ij      = ijb,ije
    173128       alpha(ij,llm) = 0.
    174129       beta (ij,llm) = 1./ unpl2k
    175       ENDDO
    176 c$OMP ENDDO NOWAIT
    177 c
    178 c     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
    179 c
    180       DO l = llm -1 , 2 , -1
    181 c
    182 c$OMP DO SCHEDULE(STATIC)
    183         DO ij = ijb, ije
    184         dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
    185         alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
    186         beta (ij,l)  =   p(ij,l  ) / dellta   
    187         ENDDO
    188 c$OMP ENDDO NOWAIT
    189 c
    190       ENDDO
     130    ENDDO
     131    !$OMP ENDDO NOWAIT
     132    !
     133    !     ... Calcul des coeff. alpha et beta  pour l = llm-1  a l = 2 ...
     134    !
     135    DO l = llm -1 , 2 , -1
     136       !
     137       !$OMP DO SCHEDULE(STATIC)
     138       DO ij = ijb, ije
     139          dellta = p(ij,l)* unpl2k + p(ij,l+1)* ( beta(ij,l+1)-unpl2k )
     140          alpha(ij,l)  = - p(ij,l+1) / dellta * alpha(ij,l+1)
     141          beta (ij,l)  =   p(ij,l  ) / dellta   
     142       ENDDO
     143       !$OMP ENDDO NOWAIT
     144    ENDDO
    191145
    192 c
    193 c  ***********************************************************************
    194 c     .....  Calcul de pk pour la couche 1 , pres du sol  ....
    195 c
    196 c$OMP DO SCHEDULE(STATIC)
    197       DO   ij   = ijb, ije
    198        pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  /
    199      *    (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
    200       ENDDO
    201 c$OMP ENDDO NOWAIT
    202 c
    203 c    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
    204 c
    205       DO l = 2, llm
    206 c$OMP DO SCHEDULE(STATIC)
    207         DO   ij   = ijb, ije
    208          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
    209         ENDDO
    210 c$OMP ENDDO NOWAIT       
    211       ENDDO
    212 c
    213 c
    214 c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    215       DO l = 1, llm
    216 c$OMP DO SCHEDULE(STATIC)
    217          DO   ij   = ijb, ije
    218            pkf(ij,l)=pk(ij,l)
    219          ENDDO
    220 c$OMP ENDDO NOWAIT             
    221       ENDDO
     146    !  ***********************************************************************
     147    !     .....  Calcul de pk pour la couche 1 , pres du sol  ....
     148    !
     149    !$OMP DO SCHEDULE(STATIC)
     150    DO   ij   = ijb, ije
     151       pk(ij,1) = ( p(ij,1)*pks(ij) - 0.5*alpha(ij,2)*p(ij,2) )  / &
     152            (  p(ij,1)* (1.+kappa) + 0.5*( beta(ij,2)-unpl2k )* p(ij,2) )
     153    ENDDO
     154    !$OMP ENDDO NOWAIT
     155    !
     156    !    ..... Calcul de pk(ij,l) , pour l = 2 a l = llm  ........
     157    !
     158    DO l = 2, llm
     159       !$OMP DO SCHEDULE(STATIC)
     160       DO   ij   = ijb, ije
     161          pk(ij,l) = alpha(ij,l) + beta(ij,l) * pk(ij,l-1)
     162       ENDDO
     163       !$OMP ENDDO NOWAIT       
     164    ENDDO
    222165
    223 c$OMP BARRIER
    224      
    225       jjb=jj_begin
    226       jje=jj_end
    227       CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
    228      
     166    if (present(pkf)) then
     167       !    calcul de pkf
    229168
    230       RETURN
    231       END
     169       DO l = 1, llm
     170          !$OMP DO SCHEDULE(STATIC)
     171          DO   ij   = ijb, ije
     172             pkf(ij,l)=pk(ij,l)
     173          ENDDO
     174          !$OMP ENDDO NOWAIT             
     175       ENDDO
     176
     177       !$OMP BARRIER
     178
     179       jjb=jj_begin
     180       jje=jj_end
     181       CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
     182    end if
     183
     184  END SUBROUTINE exner_hyb_p
     185
     186end module exner_hyb_p_m

  • LMDZ5/trunk/libf/dyn3dpar/exner_milieu_p_m.F90

    r1992 r2021  
    1 !
    2 ! $Id $
    3 !
    4       SUBROUTINE  exner_milieu_p ( ngrid, ps, p,beta, pks, pk, pkf )
    5 c
    6 c     Auteurs :  F. Forget , Y. Wanherdrick
    7 c P.Le Van  , Fr. Hourdin  .
    8 c    ..........
    9 c
    10 c    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
    11 c    .... beta, pks,pk,pkf   sont des argum.de sortie au sous-prog ...
    12 c
    13 c   ************************************************************************
    14 c    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
    15 c    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
    16 c    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
    17 c   ************************************************************************
    18 c    .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
    19 c    la pression et la fonction d'Exner  au  sol  .
    20 c
    21 c     WARNING : CECI est une version speciale de exner_hyb originale
    22 c               Utilise dans la version martienne pour pouvoir
    23 c               tourner avec des coordonnees verticales complexe
    24 c              => Il ne verifie PAS la condition la proportionalite en
    25 c              energie totale/ interne / potentielle (F.Forget 2001)
    26 c    ( voir note de Fr.Hourdin )  ,
    27 c
    28       USE parallel_lmdz
    29       IMPLICIT NONE
    30 c
    31 #include "dimensions.h"
    32 #include "paramet.h"
    33 #include "comconst.h"
    34 #include "comgeom.h"
    35 #include "comvert.h"
    36 #include "serre.h"
     1module exner_milieu_p_m
    372
    38       INTEGER  ngrid
    39       REAL p(ngrid,llmp1),pk(ngrid,llm),pkf(ngrid,llm)
    40       REAL ps(ngrid),pks(ngrid), beta(ngrid,llm)
     3  IMPLICIT NONE
    414
    42 c    .... variables locales   ...
     5contains
    436
    44       INTEGER l, ij
    45       REAL dum1
     7  SUBROUTINE  exner_milieu_p ( ngrid, ps, p, pks, pk, pkf )
     8    !
     9    !     Auteurs :  F. Forget , Y. Wanherdrick
     10    ! P.Le Van  , Fr. Hourdin  .
     11    !    ..........
     12    !
     13    !    ....  ngrid, ps,p             sont des argum.d'entree  au sous-prog ...
     14    !    ....  pks,pk,pkf   sont des argum.de sortie au sous-prog ...
     15    !
     16    !   ************************************************************************
     17    !    Calcule la fonction d'Exner pk = Cp * (p/preff) ** kappa , aux milieux des
     18    !    couches .   Pk(l) sera calcule aux milieux  des couches l ,entre les
     19    !    pressions p(l) et p(l+1) ,definis aux interfaces des llm couches .
     20    !   ************************************************************************
     21    !  .. N.B : Au sommet de l'atmosphere,  p(llm+1) = 0. , et ps et pks sont
     22    !    la pression et la fonction d'Exner  au  sol  .
     23    !
     24    !     WARNING : CECI est une version speciale de exner_hyb originale
     25    !               Utilise dans la version martienne pour pouvoir
     26    !               tourner avec des coordonnees verticales complexe
     27    !              => Il ne verifie PAS la condition la proportionalite en
     28    !              energie totale/ interne / potentielle (F.Forget 2001)
     29    !    ( voir note de Fr.Hourdin )  ,
     30    !
     31    USE parallel_lmdz
     32    !
     33    include "dimensions.h"
     34    include "paramet.h"
     35    include "comconst.h"
     36    include "comgeom.h"
     37    include "comvert.h"
     38    include "serre.h"
    4639
    47       REAL ppn(iim),pps(iim)
    48       REAL xpn, xps
    49       REAL SSUM
    50       EXTERNAL SSUM
    51       INTEGER ije,ijb,jje,jjb
    52       logical,save :: firstcall=.true.
    53 !$OMP THREADPRIVATE(firstcall)
    54       character(len=*),parameter :: modname="exner_milieu_p"
     40    INTEGER  ngrid
     41    REAL p(ngrid,llmp1),pk(ngrid,llm)
     42    REAL, optional:: pkf(ngrid,llm)
     43    REAL ps(ngrid),pks(ngrid)
    5544
    56       ! Sanity check
    57       if (firstcall) then
    58         ! sanity checks for Shallow Water case (1 vertical layer)
    59         if (llm.eq.1) then
     45    !    .... variables locales   ...
     46
     47    INTEGER l, ij
     48    REAL dum1
     49
     50    logical,save :: firstcall=.true.
     51    !$OMP THREADPRIVATE(firstcall)
     52    character(len=*),parameter :: modname="exner_milieu_p"
     53
     54    ! Sanity check
     55    if (firstcall) then
     56       ! sanity checks for Shallow Water case (1 vertical layer)
     57       if (llm.eq.1) then
    6058          if (kappa.ne.1) then
    61             call abort_gcm(modname,
    62      &      "kappa!=1 , but running in Shallow Water mode!!",42)
     59             call abort_gcm(modname, &
     60                  "kappa!=1 , but running in Shallow Water mode!!",42)
    6361          endif
    6462          if (cpp.ne.r) then
    65             call abort_gcm(modname,
    66      &      "cpp!=r , but running in Shallow Water mode!!",42)
     63             call abort_gcm(modname, &
     64                  "cpp!=r , but running in Shallow Water mode!!",42)
    6765          endif
    68         endif ! of if (llm.eq.1)
     66       endif ! of if (llm.eq.1)
    6967
    70         firstcall=.false.
    71       endif ! of if (firstcall)
    72      
    73 c$OMP BARRIER
     68       firstcall=.false.
     69    endif ! of if (firstcall)
    7470
    75 ! Specific behaviour for Shallow Water (1 vertical layer) case
    76       if (llm.eq.1) then
    77              
    78         ! Compute pks(:),pk(:),pkf(:)
    79         ijb=ij_begin
    80         ije=ij_end
    81 !$OMP DO SCHEDULE(STATIC)
    82         DO ij=ijb, ije
    83           pks(ij)=(cpp/preff)*ps(ij)
     71    !$OMP BARRIER
     72
     73    ! Specific behaviour for Shallow Water (1 vertical layer) case:
     74    if (llm.eq.1) then
     75
     76       ! Compute pks(:),pk(:),pkf(:)
     77       ijb=ij_begin
     78       ije=ij_end
     79       !$OMP DO SCHEDULE(STATIC)
     80       DO ij=ijb, ije
     81          pks(ij) = (cpp/preff) * ps(ij)
    8482          pk(ij,1) = .5*pks(ij)
    85           pkf(ij,1)=pk(ij,1)
    86         ENDDO
    87 !$OMP ENDDO
     83          if (present(pkf)) pkf(ij,1)=pk(ij,1)
     84       ENDDO
     85       !$OMP ENDDO
    8886
    89 !$OMP MASTER
    90       if (pole_nord) then
    91         DO  ij   = 1, iim
    92           ppn(ij) = aire(   ij   ) * pks(  ij     )
    93         ENDDO
    94         xpn      = SSUM(iim,ppn,1) /apoln
    95  
    96         DO ij   = 1, iip1
    97           pks(   ij     )  =  xpn
    98           pk(ij,1) = .5*pks(ij)
    99           pkf(ij,1)=pk(ij,1)
    100         ENDDO
    101       endif
    102      
    103       if (pole_sud) then
    104         DO  ij   = 1, iim
    105           pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    106         ENDDO
    107         xps      = SSUM(iim,pps,1) /apols
    108  
    109         DO ij   = 1, iip1
    110           pks( ij+ip1jm )  =  xps
    111           pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
    112           pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
    113         ENDDO
    114       endif
    115 !$OMP END MASTER
    116 !$OMP BARRIER
    117         jjb=jj_begin
    118         jje=jj_end
    119         CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
     87       !$OMP BARRIER
     88       if (present(pkf)) then
     89          jjb=jj_begin
     90          jje=jj_end
     91          CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
     92       end if
    12093
    121         ! our work is done, exit routine
    122         return
    123       endif ! of if (llm.eq.1)
     94       ! our work is done, exit routine
     95       return
     96    endif ! of if (llm.eq.1)
    12497
    125 !!!! General case:
     98    ! General case:
    12699
    127 c     -------------
    128 c     Calcul de pks
    129 c     -------------
    130    
    131       ijb=ij_begin
    132       ije=ij_end
     100    !     -------------
     101    !     Calcul de pks
     102    !     -------------
    133103
    134 c$OMP DO SCHEDULE(STATIC)
    135       DO   ij  = ijb, ije
    136         pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
    137       ENDDO
    138 c$OMP ENDDO
    139 c Synchro OPENMP ici
     104    ijb=ij_begin
     105    ije=ij_end
    140106
    141 c$OMP MASTER
    142       if (pole_nord) then
    143         DO  ij   = 1, iim
    144           ppn(ij) = aire(   ij   ) * pks(  ij     )
    145         ENDDO
    146         xpn      = SSUM(iim,ppn,1) /apoln
    147  
    148         DO ij   = 1, iip1
    149           pks(   ij     )  =  xpn
    150         ENDDO
    151       endif
    152      
    153       if (pole_sud) then
    154         DO  ij   = 1, iim
    155           pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
    156         ENDDO
    157         xps      = SSUM(iim,pps,1) /apols
    158  
    159         DO ij   = 1, iip1
    160           pks( ij+ip1jm )  =  xps
    161         ENDDO
    162       endif
    163 c$OMP END MASTER
    164 c$OMP BARRIER
    165 c
    166 c
    167 c    .... Calcul de pk  pour la couche l
    168 c    --------------------------------------------
    169 c
    170       dum1 = cpp * (2*preff)**(-kappa)
    171       DO l = 1, llm-1
    172 c$OMP DO SCHEDULE(STATIC)
    173         DO   ij   = ijb, ije
    174          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
    175         ENDDO
    176 c$OMP ENDDO NOWAIT       
    177       ENDDO
     107    !$OMP DO SCHEDULE(STATIC)
     108    DO   ij  = ijb, ije
     109       pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
     110    ENDDO
     111    !$OMP ENDDO
     112    ! Synchro OPENMP ici
    178113
    179 c    .... Calcul de pk  pour la couche l = llm ..
    180 c    (on met la meme distance (en log pression)  entre Pk(llm)
    181 c    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
     114    !$OMP BARRIER
     115    !
     116    !
     117    !    .... Calcul de pk  pour la couche l
     118    !    --------------------------------------------
     119    !
     120    dum1 = cpp * (2*preff)**(-kappa)
     121    DO l = 1, llm-1
     122       !$OMP DO SCHEDULE(STATIC)
     123       DO   ij   = ijb, ije
     124          pk(ij,l) = dum1 * (p(ij,l) + p(ij,l+1))**kappa
     125       ENDDO
     126       !$OMP ENDDO NOWAIT
     127    ENDDO
    182128
    183 c$OMP DO SCHEDULE(STATIC)
    184       DO   ij   = ijb, ije
    185          pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
    186       ENDDO
    187 c$OMP ENDDO NOWAIT       
     129    !    .... Calcul de pk  pour la couche l = llm ..
     130    !    (on met la meme distance (en log pression)  entre Pk(llm)
     131    !    et Pk(llm -1) qu'entre Pk(llm-1) et Pk(llm-2)
    188132
     133    !$OMP DO SCHEDULE(STATIC)
     134    DO   ij   = ijb, ije
     135       pk(ij,llm) = pk(ij,llm-1)**2 / pk(ij,llm-2)
     136    ENDDO
     137    !$OMP ENDDO NOWAIT       
    189138
    190 c    calcul de pkf
    191 c    -------------
    192 c      CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
    193       DO l = 1, llm
    194 c$OMP DO SCHEDULE(STATIC)
    195          DO   ij   = ijb, ije
    196            pkf(ij,l)=pk(ij,l)
    197          ENDDO
    198 c$OMP ENDDO NOWAIT             
    199       ENDDO
     139    if (present(pkf)) then
     140       !    calcul de pkf
    200141
    201 c$OMP BARRIER
    202      
    203       jjb=jj_begin
    204       jje=jj_end
    205       CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
    206      
    207 c    EST-CE UTILE ?? : calcul de beta
    208 c    --------------------------------
    209       DO l = 2, llm
    210 c$OMP DO SCHEDULE(STATIC)
    211         DO   ij   = ijb, ije
    212           beta(ij,l) = pk(ij,l) / pk(ij,l-1)   
    213         ENDDO
    214 c$OMP ENDDO NOWAIT             
    215       ENDDO
     142       DO l = 1, llm
     143          !$OMP DO SCHEDULE(STATIC)
     144          DO   ij   = ijb, ije
     145             pkf(ij,l)=pk(ij,l)
     146          ENDDO
     147          !$OMP ENDDO NOWAIT
     148       ENDDO
    216149
    217       RETURN
    218       END
     150       !$OMP BARRIER
     151
     152       jjb=jj_begin
     153       jje=jj_end
     154       CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
     155    end if
     156
     157  END SUBROUTINE exner_milieu_p
     158
     159end module exner_milieu_p_m

  • LMDZ5/trunk/libf/dyn3dpar/gcm.F

    r1939 r2021  
    9999      REAL ps(ip1jmp1)                       ! pression  au sol
    100100c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
    101 c      REAL pks(ip1jmp1)                      ! exner au  sol
    102 c      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    103 c      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    104101      REAL masse(ip1jmp1,llm)                ! masse d'air
    105102      REAL phis(ip1jmp1)                     ! geopotentiel au sol
    … …  
    125122      data call_iniphys/.true./
    126123
    127 c      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
    128124c+jld variables test conservation energie
    129125c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)

  • LMDZ5/trunk/libf/dyn3dpar/guide_p_mod.F90

    r1907 r2021  
    328328!=======================================================================
    329329  SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
     330    use exner_hyb_p_m, only: exner_hyb_p
     331    use exner_milieu_p_m, only: exner_milieu_p
    330332    USE parallel_lmdz
    331333    USE control_mod
    … …  
    349351    REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
    350352    ! Variables pour fonction Exner (P milieu couche)
    351     REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
    352     REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
     353    REAL, DIMENSION (iip1,jjp1,llm)    :: pk
    353354    REAL, DIMENSION (iip1,jjp1)        :: pks   
    354355    REAL                               :: unskap
    … …  
    493494        CALL pression_p( ip1jmp1, ap, bp, ps, p )
    494495        if (pressure_exner) then
    495           CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
     496          CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk)
    496497        else
    497           CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
     498          CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk)
    498499        endif
    499500        unskap=1./kappa
    … …  
    689690!=======================================================================
    690691  SUBROUTINE guide_interp(psi,teta)
     692    use exner_hyb_p_m, only: exner_hyb_p
     693    use exner_milieu_p_m, only: exner_milieu_p
    691694  USE parallel_lmdz
    692695  USE mod_hallo
    … …  
    713716  REAL, DIMENSION (iip1,jjm,llm)     :: pbary
    714717  ! Variables pour fonction Exner (P milieu couche)
    715   REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
    716   REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
     718  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
    717719  REAL, DIMENSION (iip1,jjp1)        :: pks   
    718720  REAL                               :: unskap
    … …  
    793795        CALL pression_p( ip1jmp1, ap, bp, psi, p )
    794796        if (pressure_exner) then
    795           CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
     797          CALL exner_hyb_p(ip1jmp1,psi,p,pks,pk)
    796798        else
    797           CALL exner_milieu_p(ip1jmp1,psi,p,beta,pks,pk,pkf)
     799          CALL exner_milieu_p(ip1jmp1,psi,p,pks,pk)
    798800        endif
    799801        unskap=1./kappa

  • LMDZ5/trunk/libf/dyn3dpar/iniacademic.F90

    r1907 r2021  
    44SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
    55
     6  use exner_hyb_m, only: exner_hyb
     7  use exner_milieu_m, only: exner_milieu
    68  USE filtreg_mod
    79  USE infotrac, ONLY : nqtot
    … …  
    5456  REAL pks(ip1jmp1)                      ! exner au  sol
    5557  REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    56   REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
    5758  REAL phi(ip1jmp1,llm)                  ! geopotentiel
    5859  REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
    … …  
    7071  integer idum
    7172
    72   REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
     73  REAL zdtvr
    7374 
    7475  character(len=*),parameter :: modname="iniacademic"
    … …  
    223224        CALL pression ( ip1jmp1, ap, bp, ps, p       )
    224225        if (pressure_exner) then
    225           CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    226         else
    227           call exner_milieu(ip1jmp1,ps,p,beta,pks,pk,pkf)
     226          CALL exner_hyb( ip1jmp1, ps, p, pks, pk )
     227        else
     228          call exner_milieu(ip1jmp1,ps,p,pks,pk)
    228229        endif
    229230        CALL massdair(p,masse)

  • LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F

    r1988 r2021  
    88     &                    time_0)
    99
     10      use exner_hyb_m, only: exner_hyb
     11      use exner_milieu_m, only: exner_milieu
     12      use exner_hyb_p_m, only: exner_hyb_p
     13      use exner_milieu_p_m, only: exner_milieu_p
    1014       USE misc_mod
    1115       USE parallel_lmdz
    … …  
    149153      character*10 string10
    150154
    151       REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
    152155      REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
    153156
    … …  
    241244      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    242245      if (pressure_exner) then
    243         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     246        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
    244247      else
    245         CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     248        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
    246249      endif
    247250c$OMP END MASTER
    … …  
    705708c$OMP BARRIER
    706709         if (pressure_exner) then
    707            CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
     710           CALL exner_hyb_p(  ip1jmp1, ps, p,pks, pk, pkf )
    708711         else
    709            CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
     712           CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
    710713         endif
    711714c$OMP BARRIER
    … …  
    918921c$OMP BARRIER
    919922          if (pressure_exner) then
    920             CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
     923            CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk,pkf)
    921924          else
    922             CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
     925            CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk,pkf)
    923926          endif
    924927c$OMP BARRIER
    … …  
    10591062c$OMP BARRIER
    10601063        if (pressure_exner) then
    1061           CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     1064          CALL exner_hyb_p( ip1jmp1, ps, p, pks, pk, pkf )
    10621065        else
    1063           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
     1066          CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
    10641067        endif
    10651068c$OMP BARRIER

  • LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90

    r1907 r2021  
    2929  USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
    3030  USE indice_sol_mod
     31  use exner_hyb_m, only: exner_hyb
     32  use exner_milieu_m, only: exner_milieu
     33  use test_disvert_m, only: test_disvert
    3134#endif
    3235  IMPLICIT NONE
    … …  
    7477  CHARACTER(LEN=80)                        :: x, fmt
    7578  INTEGER                                  :: i, j, l, ji
    76   REAL,    DIMENSION(iip1,jjp1,llm)        :: alpha, beta, pk, pls, y
     79  REAL,    DIMENSION(iip1,jjp1,llm)        :: pk, pls, y
    7780  REAL,    DIMENSION(ip1jmp1)              :: pks
    7881
    … …  
    150153
    151154  CALL iniconst()
     155  call test_disvert
    152156  CALL inigeom()
    153157
    … …  
    253257  CALL pression(ip1jmp1, ap, bp, psol, p3d)
    254258  if (pressure_exner) then
    255     CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
     259    CALL exner_hyb(ip1jmp1, psol, p3d, pks, pk)
    256260  else
    257     CALL exner_milieu(ip1jmp1,psol,p3d,beta,pks,pk,y)
     261    CALL exner_milieu(ip1jmp1,psol,p3d, pks,pk)
    258262  endif
    259263  pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)
Note: See TracChangeset for help on using the changeset viewer.