Ignore:
Timestamp:
Apr 25, 2014, 12:20:14 PM (11 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.

File:
1 moved

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.