Ignore:
Timestamp:
Jun 26, 2014, 6:07:05 PM (11 years ago)
Author:
emillour
Message:

Common dynamics:
Some updates to keep up with LMDZ5 Earth model evolution
(up to LMDZ5 rev 2070). See file "DOC/chantiers/commit_importants.log"
for detailed list of changes.
Note that the updates of exner* routines change (as expected) results
at numerical roundoff level.
EM

Location:
trunk/LMDZ.COMMON/libf/dyn3dpar
Files:
9 edited
2 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/dyn3dpar/calfis_p.F

    r1256 r1302  
    225225      save unskap
    226226
    227 cIM diagnostique PVteta, Amip2
    228       INTEGER,PARAMETER :: ntetaSTD=3
    229       REAL,SAVE :: rtetaSTD(ntetaSTD)=(/350.,380.,405./) ! Earth-specific, beware !!
    230       REAL PVteta(klon,ntetaSTD)
    231      
    232      
    233227      REAL SSUM
    234228
    … …  
    267261      klon=klon_mpi
    268262     
    269       PVteta(:,:)=0.
    270            
    271263      IF ( firstcal )  THEN
    272264        debut = .TRUE.
    … …  
    684676      endif
    685677
    686 
    687       IF (is_sequential.and.(planet_type=="earth")) THEN
    688 #ifdef CPP_EARTH
    689 ! PVtheta calls tetalevel, which is in the physics
    690 cIM calcul PV a teta=350, 380, 405K
    691         CALL PVtheta(ngridmx,llm,pucov,pvcov,pteta,
    692      $           ztfi,zplay,zplev,
    693      $           ntetaSTD,rtetaSTD,PVteta)
    694 c
    695 #endif
    696       ENDIF
    697 
    698678c On change de grille, dynamique vers physiq, pour le flux de masse verticale
    699679      CALL gr_dyn_fi_p(llm,iip1,jjp1,klon,flxw,flxwfi)
    … …  
    923903     .             zdqfi_omp,
    924904     .             zdpsrf_omp,
    925 cIM diagnostique PVteta, Amip2         
    926      .             pducov,
    927      .             PVteta)
     905     .             pducov)
    928906
    929907      else if ( planet_type=="generic" ) then

  • trunk/LMDZ.COMMON/libf/dyn3dpar/ce0l.F90

    r1019 r1302  
    11!
    2 ! $Id: ce0l.F90 1615 2012-02-10 15:42:26Z emillour $
     2! $Id: ce0l.F90 1984 2014-02-18 09:59:29Z emillour $
    33!
    44!-------------------------------------------------------------------------------
    … …  
    115115  END IF
    116116
    117   IF (grilles_gcm_netcdf) THEN
    118      WRITE(lunout,'(//)')
    119      WRITE(lunout,*) '  ***************************  '
    120      WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
    121      WRITE(lunout,*) '  ***************************  '
    122      WRITE(lunout,'(//)')
    123      CALL grilles_gcm_netcdf_sub(masque,phis)
    124   END IF
     117  WRITE(lunout,'(//)')
     118  WRITE(lunout,*) '  ***************************  '
     119  WRITE(lunout,*) '  ***  grilles_gcm_netcdf ***  '
     120  WRITE(lunout,*) '  ***************************  '
     121  WRITE(lunout,'(//)')
     122  CALL grilles_gcm_netcdf_sub(masque,phis)
    125123 
    126124#ifdef CPP_MPI
    … …  
    137135!
    138136!-------------------------------------------------------------------------------
     137

  • trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F

    r1300 r1302  
    22! $Id: conf_gcm.F 1438 2010-10-08 10:19:34Z jghattas $
    33!
    4 c
    5 c
     4!
     5!
    66      SUBROUTINE conf_gcm( tapedef, etatinit )
    7 c
     7!
    88#ifdef CPP_IOIPSL
    99      use IOIPSL
    … …  
    2020      use sponge_mod_p, only: callsponge,mode_sponge,nsponge,tetasponge
    2121      IMPLICIT NONE
    22 c-----------------------------------------------------------------------
    23 c     Auteurs :   L. Fairhead , P. Le Van  .
    24 c
    25 c     Arguments :
    26 c
    27 c     tapedef   :
    28 c     etatinit  :     = TRUE   , on ne  compare pas les valeurs des para-
    29 c     -metres  du zoom  avec  celles lues sur le fichier start .
    30 c
     22!-----------------------------------------------------------------------
     23!     Auteurs :   L. Fairhead , P. Le Van  .
     24!
     25!     Arguments :
     26!
     27!     tapedef   :
     28!     etatinit  :     = TRUE   , on ne  compare pas les valeurs des para-
     29!     -metres  du zoom  avec  celles lues sur le fichier start .
     30!
    3131       LOGICAL etatinit
    3232       INTEGER tapedef
    3333
    34 c   Declarations :
    35 c   --------------
     34!   Declarations :
     35!   --------------
    3636#include "dimensions.h"
    3737#include "paramet.h"
    … …  
    4545! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
    4646! #include "clesphys.h"
    47 c
    48 c
    49 c   local:
    50 c   ------
     47!
     48!
     49!   local:
     50!   ------
    5151
    5252      CHARACTER ch1*72,ch2*72,ch3*72,ch4*12
    … …  
    6060      integer,external :: OMP_GET_NUM_THREADS
    6161#endif     
    62 c
    63 c  -------------------------------------------------------------------
    64 c
    65 c       .........     Version  du 29/04/97       ..........
    66 c
    67 c   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
    68 c      tetatemp   ajoutes  pour la dissipation   .
    69 c
    70 c   Autre parametre ajoute en fin de liste de tapedef : ** fxyhypb **
    71 c
    72 c  Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
    73 c    Sinon , choix de fxynew  , a derivee sinusoidale  ..
    74 c
    75 c   ......  etatinit = . TRUE. si defrun  est appele dans ETAT0_LMD  ou
    76 c         LIMIT_LMD  pour l'initialisation de start.dat (dic) et
    77 c                de limit.dat ( dic)                        ...........
    78 c           Sinon  etatinit = . FALSE .
    79 c
    80 c   Donc etatinit = .F.  si on veut comparer les valeurs de  grossismx ,
    81 c    grossismy,clon,clat, fxyhypb  lues sur  le fichier  start  avec
    82 c   celles passees  par run.def ,  au debut du gcm, apres l'appel a
    83 c    lectba . 
    84 c   Ces parmetres definissant entre autres la grille et doivent etre
    85 c   pareils et coherents , sinon il y aura  divergence du gcm .
    86 c
    87 c-----------------------------------------------------------------------
    88 c   initialisations:
    89 c   ----------------
     62!
     63!  -------------------------------------------------------------------
     64!
     65!       .........     Version  du 29/04/97       ..........
     66!
     67!   Nouveaux parametres nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,
     68!      tetatemp   ajoutes  pour la dissipation   .
     69!
     70!   Autre parametre ajoute en fin de liste de tapedef : ** fxyhypb **
     71!
     72!  Si fxyhypb = .TRUE. , choix de la fonction a derivee tangente hyperb.
     73!    Sinon , choix de fxynew  , a derivee sinusoidale  ..
     74!
     75!   ......  etatinit = . TRUE. si defrun  est appele dans ETAT0_LMD  ou
     76!         LIMIT_LMD  pour l'initialisation de start.dat (dic) et
     77!                de limit.dat ( dic)                        ...........
     78!           Sinon  etatinit = . FALSE .
     79!
     80!   Donc etatinit = .F.  si on veut comparer les valeurs de  grossismx ,
     81!    grossismy,clon,clat, fxyhypb  lues sur  le fichier  start  avec
     82!   celles passees  par run.def ,  au debut du gcm, apres l'appel a
     83!    lectba . 
     84!   Ces parmetres definissant entre autres la grille et doivent etre
     85!   pareils et coherents , sinon il y aura  divergence du gcm .
     86!
     87!-----------------------------------------------------------------------
     88!   initialisations:
     89!   ----------------
    9090
    9191!Config  Key  = lunout
    … …  
    290290       CALL getin('dissip_period',dissip_period)
    291291
    292 ccc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
    293 ccc
     292!cc  ....   P. Le Van , modif le 29/04/97 .pour la dissipation  ...
     293!cc
    294294
    295295!Config  Key  = lstardis
    … …  
    456456       CALL getin('ok_guide',ok_guide)
    457457
    458 c    ...............................................................
     458!     ...............................................................
    459459
    460460!Config  Key  =  read_start
    … …  
    632632      CALL getin('ok_etat0',ok_etat0)
    633633
    634 !Config  Key  = grilles_gcm_netcdf
    635 !Config  Desc = creation de fichier grilles_gcm.nc dans create_etat0_limit
    636 !Config  Def  = n
    637       grilles_gcm_netcdf = .FALSE.
    638       CALL getin('grilles_gcm_netcdf',grilles_gcm_netcdf)
    639 
    640 c----------------------------------------
    641 c Parameters for zonal averages in the case of Titan
     634!----------------------------------------
     635! Parameters for zonal averages in the case of Titan
    642636      moyzon_mu = .false.
    643637      moyzon_ch = .false.
    … …  
    646640       CALL getin('moyzon_ch', moyzon_ch)
    647641      endif
    648 c----------------------------------------
    649 
    650 c----------------------------------------
    651 ccc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
    652 c     .........   (  modif  le 17/04/96 )   .........
    653 c
    654 C ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
    655 c
    656 c----------------------------------------
     642!----------------------------------------
     643
     644!----------------------------------------
     645!cc  ....   P. Le Van , ajout  le 7/03/95 .pour le zoom ...
     646!     .........   (  modif  le 17/04/96 )   .........
     647!
     648! ZOOM PARAMETERS ... the ones read in start.nc prevail anyway ! (SL, 2012)
     649!
     650!----------------------------------------
    657651      IF( etatinit ) then
    658652
    … …  
    707701
    708702      write(lunout,*)'conf_gcm: alphax alphay ',alphax,alphay
    709 c
    710 c    alphax et alphay sont les anciennes formulat. des grossissements
    711 c
    712 c
     703!
     704!    alphax et alphay sont les anciennes formulat. des grossissements
     705!
     706!
    713707
    714708!Config  Key  = fxyhypb
    … …  
    758752       ysinus = .TRUE.
    759753       CALL getin('ysinus',ysinus)
    760 c
    761 c----------------------------------------
     754!
     755!----------------------------------------
    762756       else ! etatinit=false
    763 c----------------------------------------
     757!----------------------------------------
    764758
    765759!Config  Key  = clon
    … …  
    779773       CALL getin('clat',clatt)
    780774
    781 c
    782 c
     775!
     776!
    783777      IF( ABS(clat - clatt).GE. 0.001 )  THEN
    784778        write(lunout,*)'conf_gcm: La valeur de clat passee par run.def',
    … …  
    834828
    835829      write(lunout,*)'conf_gcm: alphax alphay',alphax,alphay
    836 c
    837 c    alphax et alphay sont les anciennes formulat. des grossissements
    838 c
    839 c
     830!
     831!    alphax et alphay sont les anciennes formulat. des grossissements
     832!
     833!
    840834
    841835!Config  Key  = fxyhypb
    … …  
    862856         ENDIF
    863857      ENDIF
    864 c
     858!
    865859!Config  Key  = dzoomx
    866860!Config  Desc = extension en longitude
    … …  
    925919      ENDIF
    926920
    927 cc
     921!c
    928922      IF( .NOT.fxyhypb  )  THEN
    929923
    … …  
    955949
    956950      endif ! etatinit
    957 c----------------------------------------
     951!----------------------------------------
    958952
    959953
    … …  
    10091003      write(lunout,*)' ok_limit = ', ok_limit
    10101004      write(lunout,*)' ok_etat0 = ', ok_etat0
    1011       write(lunout,*)' grilles_gcm_netcdf = ', grilles_gcm_netcdf
    10121005      if (planet_type=="titan") then
    10131006       write(lunout,*)' moyzon_mu = ', moyzon_mu

  • trunk/LMDZ.COMMON/libf/dyn3dpar/exner_hyb_p_m.F90

    r1299 r1302  
    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
     187

  • trunk/LMDZ.COMMON/libf/dyn3dpar/exner_milieu_p_m.F90

    r1299 r1302  
    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,ijb,ije,jjb,jje
     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
     160

  • trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F

    r1300 r1302  
    516516
    517517
    518       day_end = day_ini + nday
     518      if (nday>=0) then ! standard case
     519        day_end=day_ini+nday
     520      else ! special case when nday <0, run -nday dynamical steps
     521        day_end=day_ini-nday/day_step
     522      endif
    519523      if (less1day) then
    520524        day_end=day_ini+floor(time_0+fractday)

  • trunk/LMDZ.COMMON/libf/dyn3dpar/guide_p_mod.F90

    r1300 r1302  
    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
    … …  
    491492    IF (f_out) THEN
    492493!       Calcul niveaux pression milieu de couches
    493         CALL pression_p( ip1jmp1, ap, bp, ps, p )
    494         if (pressure_exner) then
    495           CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
    496         else
    497           CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
     494        CALL pression_p( ip1jmp1, ap, bp, ps, p )
     495        if (pressure_exner) then
     496          CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk)
     497        else
     498          CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk)
    498499        endif
    499500        unskap=1./kappa
    500         DO l = 1, llm
    501             DO j=jjb_u,jje_u
    502                 DO i =1, iip1
    503                     p(i+(j-1)*iip1,l) = preff * ( pk(i,j,l)/cpp) ** unskap
    504                 ENDDO
    505             ENDDO
    506         ENDDO
    507         CALL guide_out("P",jjp1,llm,p(1:ip1jmp1,1:llm),1.)
     501        DO l = 1, llm
     502           DO j=jjb_u,jje_u
     503              DO i =1, iip1
     504                 p(i+(j-1)*iip1,l) = preff * ( pk(i,j,l)/cpp) ** unskap
     505              ENDDO
     506           ENDDO
     507        ENDDO
     508        CALL guide_out("SP",jjp1,llm,p(1:ip1jmp1,1:llm),1.)
    508509    ENDIF
    509510   
    … …  
    517518        if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
    518519        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
    519         IF (f_out) CALL guide_out("U",jjp1,llm,f_add(:,:),factt)
     520        IF (f_out) CALL guide_out("u",jjp1,llm,ucov,factt)
     521        IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)*ugui1(:,:)+tau*ugui2(:,:),factt)
     522        IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add(:,:)/factt,factt)
    520523        ucov(ijb_u:ije_u,:)=ucov(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
    521524    endif
    … …  
    529532        if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
    530533        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
    531         IF (f_out) CALL guide_out("T",jjp1,llm,f_add(:,:),factt)
     534        IF (f_out) CALL guide_out("teta",jjp1,llm,f_add(:,:)/factt,factt)
    532535        teta(ijb_u:ije_u,:)=teta(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
    533536    endif
    … …  
    541544        if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
    542545        CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
    543         IF (f_out) CALL guide_out("SP",jjp1,1,f_add(1:ip1jmp1,1),factt)
     546        IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt,factt)
    544547        ps(ijb_u:ije_u)=ps(ijb_u:ije_u)+f_add(ijb_u:ije_u,1)
    545548        CALL pression_p(ip1jmp1,ap,bp,ps,p)
    … …  
    555558        if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
    556559        CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
    557         IF (f_out) CALL guide_out("Q",jjp1,llm,f_add(:,:),factt)
     560        IF (f_out) CALL guide_out("q",jjp1,llm,f_add(:,:)/factt,factt)
    558561        q(ijb_u:ije_u,:)=q(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
    559562    endif
    … …  
    568571        if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
    569572        CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
    570         IF (f_out) CALL guide_out("V",jjm,llm,f_add(1:ip1jm,:),factt)
     573        IF (f_out) CALL guide_out("v",jjm,llm,vcov(1:ip1jm,:),factt)
     574        IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)*vgui1(:,:)+tau*vgui2(:,:),factt)
     575        IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt,factt)
    571576        vcov(ijb_v:ije_v,:)=vcov(ijb_v:ije_v,:)+f_add(ijb_v:ije_v,:)
    572577    endif
    … …  
    689694!=======================================================================
    690695  SUBROUTINE guide_interp(psi,teta)
     696    use exner_hyb_p_m, only: exner_hyb_p
     697    use exner_milieu_p_m, only: exner_milieu_p
    691698  USE parallel_lmdz
    692699  USE mod_hallo
    … …  
    713720  REAL, DIMENSION (iip1,jjm,llm)     :: pbary
    714721  ! Variables pour fonction Exner (P milieu couche)
    715   REAL, DIMENSION (iip1,jjp1,llm)    :: pk, pkf
    716   REAL, DIMENSION (iip1,jjp1,llm)    :: alpha, beta
     722  REAL, DIMENSION (iip1,jjp1,llm)    :: pk
    717723  REAL, DIMENSION (iip1,jjp1)        :: pks   
    718724  REAL                               :: unskap
    … …  
    784790    IF (guide_plevs.EQ.1) THEN
    785791        DO l=1,llm
    786             DO j=jjb_u,jje_u
    787                 DO i =1, iip1
     792            DO j=jjb_u,jje_u
     793                DO i =1, iip1
    788794                    pls(i,j,l)=(ap(l)+ap(l+1))/2.+psi(i,j)*(bp(l)+bp(l+1))/2.
    789                 ENDDO
    790             ENDDO
     795                ENDDO
     796            ENDDO
    791797        ENDDO
    792798    ELSE
    793         CALL pression_p( ip1jmp1, ap, bp, psi, p )
    794         if (pressure_exner) then
    795           CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
     799        CALL pression_p( ip1jmp1, ap, bp, psi, p )
     800        if (pressure_exner) then
     801          CALL exner_hyb_p(ip1jmp1,psi,p,pks,pk)
    796802        else
    797           CALL exner_milieu_p(ip1jmp1,psi,p,beta,pks,pk,pkf)
     803          CALL exner_milieu_p(ip1jmp1,psi,p,pks,pk)
    798804        endif
    799         unskap=1./kappa
    800         DO l = 1, llm
    801             DO j=jjb_u,jje_u
    802                 DO i =1, iip1
    803                     pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
    804                 ENDDO
    805             ENDDO
    806         ENDDO
     805        unskap=1./kappa
     806        DO l = 1, llm
     807            DO j=jjb_u,jje_u
     808                DO i =1, iip1
     809                    pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
     810                ENDDO
     811            ENDDO
     812        ENDDO
    807813    ENDIF
    808814
    … …  
    10241030        ! Calcul des nouvelles valeurs des niveaux de pression du guidage
    10251031        IF (guide_plevs.EQ.1) THEN
    1026         CALL Register_SwapFieldHallo(psnat1,psnat1,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
    1027         CALL SendRequest(Req)
    1028         CALL WaitRequest(Req)
    1029         CALL Register_SwapFieldHallo(psnat2,psnat2,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
    1030         CALL SendRequest(Req)
    1031         CALL WaitRequest(Req)
     1032        CALL Register_SwapFieldHallo(psnat1,psnat1,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
     1033        CALL SendRequest(Req)
     1034        CALL WaitRequest(Req)
     1035        CALL Register_SwapFieldHallo(psnat2,psnat2,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
     1036        CALL SendRequest(Req)
     1037        CALL WaitRequest(Req)
    10321038            DO l=1,nlevnc
    10331039                DO j=jjb_v,jje_v
    … …  
    10411047            ENDDO
    10421048        ELSE IF (guide_plevs.EQ.2) THEN
    1043         CALL Register_SwapFieldHallo(pnat1,pnat1,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
    1044         CALL SendRequest(Req)
    1045         CALL WaitRequest(Req)
    1046         CALL Register_SwapFieldHallo(pnat2,pnat2,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
    1047         CALL SendRequest(Req)
    1048         CALL WaitRequest(Req)
     1049        CALL Register_SwapFieldHallo(pnat1,pnat1,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
     1050        CALL SendRequest(Req)
     1051        CALL WaitRequest(Req)
     1052        CALL Register_SwapFieldHallo(pnat2,pnat2,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
     1053        CALL SendRequest(Req)
     1054        CALL WaitRequest(Req)
    10491055            DO l=1,nlevnc
    10501056                DO j=jjb_v,jje_v
    … …  
    18071813   
    18081814    ! Variables entree
    1809     CHARACTER, INTENT(IN)                          :: varname
     1815    CHARACTER*(*), INTENT(IN)                          :: varname
    18101816    INTEGER,   INTENT (IN)                         :: hsize,vsize
    18111817    REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
    … …  
    18171823    INTEGER       :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
    18181824    INTEGER       :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
     1825    INTEGER       :: vid_au,vid_av
     1826    INTEGER       :: l
    18191827    INTEGER, DIMENSION (3) :: dim3
    18201828    INTEGER, DIMENSION (4) :: dim4,count,start
    18211829    INTEGER                :: ierr, varid
     1830    REAL, DIMENSION (iip1,hsize,vsize) :: field2
    18221831   
    18231832    CALL gather_field(field,iip1*hsize,vsize,0)
    … …  
    18341843! Definition des dimensions
    18351844        ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu)
     1845        print*,'id_lonu 1 ',id_lonu
    18361846        ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv)
    18371847        ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu)
    … …  
    18421852! Creation des variables dimensions
    18431853        ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
     1854        print*,'id_lonu 2 ',id_lonu
    18441855        ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
    18451856        ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
    … …  
    18481859        ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
    18491860        ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
     1861        ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
     1862        ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
    18501863       
    18511864        ierr=NF_ENDDEF(nid)
    … …  
    18531866! Enregistrement des variables dimensions
    18541867#ifdef NC_DOUBLE
     1868        print*,'id_lonu DOUBLE ',id_lonu,rlonu*180./pi
    18551869        ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
    18561870        ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
    … …  
    18601874        ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
    18611875        ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
     1876        ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
     1877        ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
    18621878#else
     1879        print*,'id_lonu 3 ',id_lonu,rlonu*180./pi
    18631880        ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
    18641881        ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
    … …  
    18681885        ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
    18691886        ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
     1887        ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
     1888        ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
    18701889#endif
    18711890! --------------------------------------------------------------------
    … …  
    18751894! Pressure (GCM)
    18761895        dim4=(/id_lonv,id_latu,id_lev,id_tim/)
    1877         ierr = NF_DEF_VAR(nid,"P",NF_FLOAT,4,dim4,varid)
     1896        ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,4,dim4,varid)
    18781897! Surface pressure (guidage)
    18791898        IF (guide_P) THEN
    … …  
    18831902! Zonal wind
    18841903        IF (guide_u) THEN
     1904        print*,'id_lonu 4 ',id_lonu,varname
    18851905            dim4=(/id_lonu,id_latu,id_lev,id_tim/)
     1906            ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
     1907            ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
    18861908            ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
    18871909        ENDIF
    … …  
    18891911        IF (guide_v) THEN
    18901912            dim4=(/id_lonv,id_latv,id_lev,id_tim/)
     1913            ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
     1914            ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
    18911915            ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
    18921916        ENDIF
    … …  
    19121936    ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
    19131937
     1938    IF (varname=="SP") timestep=timestep+1
     1939
     1940    IF (varname=="SP") THEN
     1941      print*,'varname=SP=',varname
     1942    ELSE
     1943      print*,'varname diff SP=',varname
     1944    ENDIF
     1945
     1946
     1947    ierr = NF_INQ_VARID(nid,varname,varid)
    19141948    SELECT CASE (varname)
    1915     CASE ("P")
    1916         timestep=timestep+1
    1917         ierr = NF_INQ_VARID(nid,"P",varid)
     1949    CASE ("SP","ps")
    19181950        start=(/1,1,1,timestep/)
    19191951        count=(/iip1,jjp1,llm,1/)
    1920 #ifdef NC_DOUBLE
    1921         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
    1922 #else
    1923         ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
    1924 #endif
    1925     CASE ("SP")
    1926         ierr = NF_INQ_VARID(nid,"ps",varid)
    1927         start=(/1,1,timestep,0/)
    1928         count=(/iip1,jjp1,1,0/)
    1929 #ifdef NC_DOUBLE
    1930         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
    1931 #else
    1932         ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
    1933 #endif
    1934     CASE ("U")
    1935         ierr = NF_INQ_VARID(nid,"ucov",varid)
     1952    CASE ("v","va","vcov")
     1953        start=(/1,1,1,timestep/)
     1954        count=(/iip1,jjm,llm,1/)
     1955    CASE DEFAULT
    19361956        start=(/1,1,1,timestep/)
    19371957        count=(/iip1,jjp1,llm,1/)
     1958    END SELECT
     1959
     1960    SELECT CASE (varname)
     1961    CASE("u","ua")
     1962        DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
     1963        field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
     1964    CASE("v","va")
     1965        DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
     1966    CASE DEFAULT
     1967        field2=field
     1968    END SELECT
     1969
    19381970#ifdef NC_DOUBLE
    1939         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
     1971    ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
    19401972#else
    1941         ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
     1973    ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field2)
    19421974#endif
    1943     CASE ("V")
    1944         ierr = NF_INQ_VARID(nid,"vcov",varid)
    1945         start=(/1,1,1,timestep/)
    1946         count=(/iip1,jjm,llm,1/)
    1947 #ifdef NC_DOUBLE
    1948         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
    1949 #else
    1950         ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
    1951 #endif
    1952     CASE ("T")
    1953         ierr = NF_INQ_VARID(nid,"teta",varid)
    1954         start=(/1,1,1,timestep/)
    1955         count=(/iip1,jjp1,llm,1/)
    1956 #ifdef NC_DOUBLE
    1957         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
    1958 #else
    1959         ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
    1960 #endif
    1961     CASE ("Q")
    1962         ierr = NF_INQ_VARID(nid,"q",varid)
    1963         start=(/1,1,1,timestep/)
    1964         count=(/iip1,jjp1,llm,1/)
    1965 #ifdef NC_DOUBLE
    1966         ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field/factt)
    1967 #else
    1968         ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field/factt)
    1969 #endif
    1970     END SELECT
    19711975 
    19721976    ierr = NF_CLOSE(nid)

  • trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F

    r1300 r1302  
    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
    … …  
    1721       USE vampir
    1822       USE timer_filtre, ONLY : print_filtre_timer
    19        USE infotrac
     23       USE infotrac, ONLY: nqtot
    2024       USE guide_p_mod, ONLY : guide_main
    2125       USE getparam
    … …  
    165169      character*10 string10
    166170
    167       REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
    168171      REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
    169172
    … …  
    236239      lafin=.false.
    237240     
    238       itaufin   = nday*day_step
     241      if (nday>=0) then
     242         itaufin   = nday*day_step
     243      else
     244         ! to run a given (-nday) number of dynamical steps
     245         itaufin   = -nday
     246      endif
    239247      if (less1day) then
    240248c MODIF VENUS: to run less than one day:
    … …  
    292300      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    293301      if (pressure_exner) then
    294         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     302        CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
    295303      else
    296         CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
     304        CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
    297305      endif
    298306c$OMP END MASTER
    … …  
    830838c$OMP BARRIER
    831839         if (pressure_exner) then
    832            CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
     840           CALL exner_hyb_p(  ip1jmp1, ps, p,pks, pk, pkf )
    833841         else
    834            CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
     842           CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
    835843         endif
     844! Compute geopotential (physics might need it)
     845         CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
    836846c$OMP BARRIER
    837847           jD_cur = jD_ref + day_ini - day_ref
    … …  
    10401050c$OMP BARRIER
    10411051          if (pressure_exner) then
    1042             CALL exner_hyb_p(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
     1052            CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk,pkf)
    10431053          else
    1044             CALL exner_milieu_p(ip1jmp1,ps,p,beta,pks,pk,pkf)
     1054            CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk,pkf)
    10451055          endif
    10461056c$OMP BARRIER
    … …  
    12021212c$OMP BARRIER
    12031213        if (pressure_exner) then
    1204           CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     1214          CALL exner_hyb_p( ip1jmp1, ps, p, pks, pk, pkf )
    12051215        else
    1206           CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
     1216          CALL exner_milieu_p( ip1jmp1, ps, p, pks, pk, pkf )
    12071217        endif
    12081218c$OMP BARRIER

  • trunk/LMDZ.COMMON/libf/dyn3dpar/logic.h

    r1056 r1302  
    1111     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
    1212     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
    13      &  ,ok_limit,ok_etat0,grilles_gcm_netcdf,hybrid                    &
     13     &  ,ok_limit,ok_etat0,hybrid                                       &
    1414     &  ,moyzon_mu,moyzon_ch
    1515
    … …  
    1919     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
    2020     &  ,read_start,ok_guide,ok_strato,tidal,ok_gradsfile               &
    21      &  ,ok_limit,ok_etat0,grilles_gcm_netcdf
     21     &  ,ok_limit,ok_etat0
    2222      logical hybrid ! vertical coordinate is hybrid if true (sigma otherwise)
    2323                     ! (only used if disvert_type==2)

  • trunk/LMDZ.COMMON/libf/dyn3dpar/mod_const_mpi.F90

    r1300 r1302  
    11!
    2 ! $Id: mod_const_mpi.F90 1700 2012-12-20 14:43:19Z lguez $
     2! $Id: mod_const_mpi.F90 2055 2014-06-04 12:33:27Z acaubel $
    33!
    44MODULE mod_const_mpi
    … …  
    1717    USE ioipsl_getincom, only: getin
    1818#endif
    19 
     19! Use of Oasis-MCT coupler
     20#ifdef CPP_OMCT
     21    USE mod_prism
     22#endif
     23#ifdef CPP_XIOS
     24    USE wxios, only: wxios_init
     25#endif
    2026    IMPLICIT NONE
    2127#ifdef CPP_MPI
    … …  
    3844#ifdef CPP_COUPLE
    3945!$OMP MASTER
    40        CALL prism_init_comp_proto (comp_id, 'lmdz.x', ierr)
     46#ifdef CPP_XIOS
     47      CALL wxios_init("LMDZ", outcom=COMM_LMDZ, type_ocean=type_ocean)
     48#else
     49       CALL prism_init_comp_proto (comp_id, 'LMDZ', ierr)
    4150       CALL prism_get_localcomm_proto(COMM_LMDZ,ierr)
     51#endif
    4252!$OMP END MASTER
    4353#endif

  • trunk/LMDZ.COMMON/libf/dyn3dpar/parallel_lmdz.F90

    r1300 r1302  
    225225#endif
    226226#ifdef CPP_COUPLE
     227! Use of Oasis-MCT coupler
     228#if defined CPP_OMCT
     229    use mod_prism
     230#else
    227231    use mod_prism_proto
    228232#endif
    229 #ifdef CPP_EARTH
    230233! Ehouarn: surface_data module is in 'phylmd' ...
    231234      use surface_data, only : type_ocean
    … …  
    252255
    253256      if (type_ocean == 'couple') then
     257#ifdef CPP_XIOS
     258    !Fermeture propre de XIOS
     259      CALL wxios_close()
     260#else
    254261#ifdef CPP_COUPLE
    255262         call prism_terminate_proto(ierr)
    … …  
    258265         endif
    259266#endif
     267#endif
    260268      else
    261269#ifdef CPP_XIOS
Note: See TracChangeset for help on using the changeset viewer.