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.

File:
1 moved

Legend:

Unmodified
Added
Removed

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