Changeset 124 for trunk


Ignore:
Timestamp:
May 19, 2011, 5:05:39 PM (14 years ago)
Author:
emillour
Message:

EM: suite mise au point discretisation verticale et quelques corrections de bugs dans le version de reference parallele.

Location:
trunk
Files:
19 edited

Legend:

Unmodified
Added
Removed

  • trunk/chantiers/commit_importants.log

    r119 r124  
    881881correction d'un certain nombre de bugs.
    882882
     883*********************
     884**** commit_v123 ****
     885*********************
     886
     887Ehouarn: suite de l'implémentation de la discrétisation verticale
     888-> insertion flag 'hybrid' pour discretisation verticale en sigma/hybrid
     889   dans disvert_noterre.F
     890-> ajout parametre 'scaleheight' (lu dans z2sig.def) dans comvert.h
     891-> quelques corrections de bug (en //), mais il faudra plus de travail
     892   pour atteindre la concordance entre versions séq. et //.
     893   

  • trunk/chantiers/meschantiers-Ehouarn.txt

    r8 r124  
    11>>> choses à faire:
    22
    3 - En priorité: faire qu'on puisse compiler un/des cas test (sans physique) pour pouvoir tester les modifs au fur et à mesure.
     3- En priorité: faire qu'on puisse compiler un/des cas test (sans physique) pour pouvoir tester les modifs au fur et à mesure. OK, c'est fait, on peut compiler/tourner sans physique
    44
    55- Uniformiser les mises à jour dyn séq. et // (commencé avec la rev. 8)
    … …  
    77- Peut-être revoir l'interface dynamique/physique ?
    88
    9 - Attention à la discrétisation verticale ...
     9- Attention à la discrétisation verticale ... 
    1010

  • trunk/documentation/disvert.tex

    r109 r124  
    3737
    3838\vspace{1cm}
    39 S\'ebastien Lebonnois
     39S\'ebastien Lebonnois, Ehouarn Millour
    4040
    4141\vspace{1cm}
    … …  
    7070This is done only once, called at the beginning from \textsf{iniconst.F}.
    7171
     72In the Earth version the vertical coordinates are hybrid (sigma-pressure),
     73and generated automaticaly (or generated from parameters read from file
     74\textsf{sigma.def}, if that file is present in the directory where the
     75gcm is run).
     76
     77In the planetary version, the vertical coordinates can be hybrid (default
     78behavior) or sigma (set using parameter "hybrid" in \textsf{run.def}; true
     79implies hybrid coordinate, false implies sigma coordinate). the distribution
     80of model levels is set from file \textsf{esasig.def} or \textsf{z2sig.def},
     81depending on which is present (in the directory where the gcm is run).
     82The first line of the  \textsf{z2sig.def} file should give the value of the
     83reference atmospheric scale height (in km), followed by the (rough estimate)
     84of the altitude (in km) of the atmospheric level (one per line of the file).
     85
    7286\item Interface pressures:
    7387computed in \textsf{caldyn0.F, caldyn.F, integrd.F, leapfrog.F}

  • trunk/libf/dyn3d/comvert.h

    r109 r124  
    77      COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),      &
    88     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
    9      &               aps(llm),bps(llm)
     9     &               aps(llm),bps(llm),scaleheight
    1010
    1111      REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig,aps,bps
     12      real scaleheight ! atmospheric (reference) scale height (km)
    1213
    1314 !-----------------------------------------------------------------------

  • trunk/libf/dyn3d/disvert_noterre.F

    r110 r124  
    44c    Nouvelle version 100% Mars !!
    55c    On l'utilise aussi pour Venus et Titan, legerment modifiee.
     6
     7#ifdef CPP_IOIPSL
     8      use IOIPSL
     9#else
     10! if not using IOIPSL, we still need to use (a local version of) getin
     11      use ioipsl_getincom
     12#endif
    613
    714      IMPLICIT NONE
    … …  
    1219#include "comconst.h"
    1320#include "logic.h"
     21#include "iniprint.h"
    1422c
    1523c=======================================================================
    … …  
    2432      INTEGER l,ll
    2533      REAL snorm
    26       REAL alpha,beta,gama,delta,deltaz,h,quoi,quand
     34      REAL alpha,beta,gama,delta,deltaz
     35      real quoi,quand
    2736      REAL zsig(llm),sig(llm+1)
    2837      INTEGER np,ierr
    … …  
    4453c-----------------------------------------------------------------------
    4554c
     55! Initializations:
    4656      pi=2.*ASIN(1.)
     57     
     58      hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates)
     59      CALL getin('hybrid',hybrid)
     60      write(lunout,*)'disvert_noterre: hybrid=',hybrid
    4761
    4862! Ouverture possible de fichiers typiquement E.T.
    … …  
    6781c        <-> energie cinetique, d'apres la note de Frederic Hourdin...
    6882
    69          PRINT*,'*****************************'
    70          PRINT*,'WARNING Lecture de esasig.def'
    71          PRINT*,'*****************************'
    72          READ(99,*) h
     83         write(lunout,*)'*****************************'
     84         write(lunout,*)'WARNING reading esasig.def'
     85         write(lunout,*)'*****************************'
     86         READ(99,*) scaleheight
    7387         READ(99,*) dz0
    7488         READ(99,*) dz1
    … …  
    7690         CLOSE(99)
    7791
    78          dz0=dz0/h
    79          dz1=dz1/h
     92         dz0=dz0/scaleheight
     93         dz1=dz1/scaleheight
    8094
    8195         sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut)
    … …  
    116130
    117131      ELSE IF(ierr4.eq.0) then
    118          PRINT*,'****************************'
    119          PRINT*,'Lecture de z2sig.def'
    120          PRINT*,'****************************'
    121 
    122          READ(99,*) h
     132         write(lunout,*)'****************************'
     133         write(lunout,*)'Reading z2sig.def'
     134         write(lunout,*)'****************************'
     135
     136         READ(99,*) scaleheight
    123137         do l=1,llm
    124138            read(99,*) zsig(l)
    … …  
    128142         sig(1) =1
    129143         do l=2,llm
    130            sig(l) = 0.5 * ( exp(-zsig(l)/h) + exp(-zsig(l-1)/h) )
     144           sig(l) = 0.5 * ( exp(-zsig(l)/scaleheight) +
     145     &                      exp(-zsig(l-1)/scaleheight) )
    131146         end do
    132147         sig(llm+1) =0
    … …  
    134149c-----------------------------------------------------------------------
    135150      ELSE
    136          write(*,*) 'didn t you forget something ??? '
    137          write(*,*) 'We need the file  z2sig.def ! (OR esasig.def) '
     151         write(lunout,*) 'didn t you forget something ??? '
     152         write(lunout,*) 'We need file  z2sig.def ! (OR esasig.def)'
    138153         stop
    139154      ENDIF
    … …  
    156171c-----------------------------------------------------------------------
    157172c
    158 c se trouvait dans Titan... Verifier...
    159 c     h = 40.
    160 
    161       if (1.eq.1) then  ! toujours hybrides
    162          write(*,*) "*******************************"
    163          write(*,*) "Systeme en coordonnees hybrides"
    164          write(*,*)
     173
     174      if (hybrid) then  ! use hybrid coordinates
     175         write(lunout,*) "*********************************"
     176         write(lunout,*) "Using hybrid vertical coordinates"
     177         write(lunout,*)
    165178c        Coordonnees hybrides avec mod
    166179         DO l = 1, llm
    … …  
    172185         bp(llmp1) = 0.
    173186         ap(llmp1) = 0.
    174       else
    175          write(*,*) "****************************"
    176          write(*,*) "Systeme en coordonnees sigma"
    177          write(*,*)
     187      else ! use sigma coordinates
     188         write(lunout,*) "********************************"
     189         write(lunout,*) "Using sigma vertical coordinates"
     190         write(lunout,*)
    178191c        Pour ne pas passer en coordonnees hybrides
    179192         DO l = 1, llm
    … …  
    186199      bp(llmp1) =   0.
    187200
    188       PRINT *,' BP '
    189       PRINT *,  bp
    190       PRINT *,' AP '
    191       PRINT *,  ap
     201      write(lunout,*)' BP '
     202      write(lunout,*)  bp
     203      write(lunout,*)' AP '
     204      write(lunout,*)  ap
    192205
    193206c     Calcul au milieu des couches :
    … …  
    197210c     (on met la meme distance (en log pression)  entre P(llm)
    198211c     et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
    199 c     Specifique.  Ce choix est spécifié ici ET dans exner_hyb.F
     212c     Specifique.  Ce choix est spécifié ici ET dans exner_milieu.F
    200213
    201214      DO l = 1, llm-1
    … …  
    204217      ENDDO
    205218     
    206 c     if (hybrid) then
     219      if (hybrid) then
    207220         aps(llm) = aps(llm-1)**2 / aps(llm-2)
    208221         bps(llm) = 0.5*(bp(llm) + bp(llm+1))
    209 c     else
    210 c        bps(llm) = bps(llm-1)**2 / bps(llm-2)
    211 c        aps(llm) = 0.
    212 c     end if
    213 
    214       PRINT *,' BPs '
    215       PRINT *,  bps
    216       PRINT *,' APs'
    217       PRINT *,  aps
     222      else
     223         bps(llm) = bps(llm-1)**2 / bps(llm-2)
     224         aps(llm) = 0.
     225      end if
     226
     227      write(lunout,*)' BPs '
     228      write(lunout,*)  bps
     229      write(lunout,*)' APs'
     230      write(lunout,*)  aps
    218231
    219232      DO l = 1, llm
    220233       presnivs(l) = aps(l)+bps(l)*preff
    221        pseudoalt(l) = -h*log(presnivs(l)/preff)
     234       pseudoalt(l) = -scaleheight*log(presnivs(l)/preff)
    222235      ENDDO
    223236
    224       PRINT *,' PRESNIVS'
    225       PRINT *,presnivs
    226       PRINT *,'Pseudo altitude des Presnivs : (avec H = ',h,' km)'
    227       PRINT *,pseudoalt
     237      write(lunout,*)' PRESNIVS'
     238      write(lunout,*)presnivs
     239      write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ',
     240     &                'height of ',scaleheight,' km)'
     241      write(lunout,*)pseudoalt
    228242
    229243c     --------------------------------------------------
    … …  
    232246c     --------------------------------------------------
    233247c     open (53,file='testhybrid.tab')
    234 c     h=15.5
     248c     scaleheight=15.5
    235249c     do iz=0,34
    236250c       z = -5 + min(iz,34-iz)
    237251c     approximation of scale height for Venus
    238 c       h = 15.5 - z/55.*10.
    239 c       ps = preff*exp(-z/h)
    240 c       zsig(1)= -h*log((aps(1) + bps(1)*ps)/preff)
     252c       scaleheight = 15.5 - z/55.*10.
     253c       ps = preff*exp(-z/scaleheight)
     254c       zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff)
    241255c       do l=2,llm
    242256c     approximation of scale height for Venus
    243257c          if (zsig(l-1).le.55.) then
    244 c             h = 15.5 - zsig(l-1)/55.*10.
     258c             scaleheight = 15.5 - zsig(l-1)/55.*10.
    245259c          else
    246 c             h = 5.5 - (zsig(l-1)-55.)/35.*2.
     260c             scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2.
    247261c          endif
    248 c          zsig(l)= zsig(l-1)-h*
     262c          zsig(l)= zsig(l-1)-scaleheight*
    249263c    .    log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps))
    250264c       end do

  • trunk/libf/dyn3d/exner_milieu.F

    r109 r124  
     1!
     2! $Id $
     3!
    14      SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
    25c
    … …  
    4447      REAL xpn, xps
    4548      REAL SSUM
    46       EXTERNAL filtreg, SSUM
     49      EXTERNAL SSUM
     50
     51      if (llm.eq.1) then
     52        ! Specific behaviour for Shallow Water (1 vertical layer) case
     53     
     54        ! Sanity checks
     55        if (kappa.ne.1) then
     56          call abort_gcm("exner_hyb",
     57     &    "kappa!=1 , but running in Shallow Water mode!!",42)
     58        endif
     59        if (cpp.ne.r) then
     60        call abort_gcm("exner_hyb",
     61     &    "cpp!=r , but running in Shallow Water mode!!",42)
     62        endif
     63       
     64        ! Compute pks(:),pk(:),pkf(:)
     65       
     66        DO   ij  = 1, ngrid
     67          pks(ij) = (cpp/preff) * ps(ij)
     68          pk(ij,1) = .5*pks(ij)
     69        ENDDO
     70       
     71        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
     72        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
     73       
     74        ! our work is done, exit routine
     75        return
     76      endif ! of if (llm.eq.1)
     77
    4778     
    4879c     -------------

  • trunk/libf/dyn3d/iniacademic.F90

    r66 r124  
    207207        enddo
    208208     enddo
    209      write(lunout,*) 'Iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
     209!     write(lunout,*) 'Iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
    210210
    211211     ! 3. Initialize fields (if necessary)
    … …  
    217217
    218218        CALL pression ( ip1jmp1, ap, bp, ps, p       )
    219         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     219        if (planet_type=="earth") then
     220          CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
     221        else
     222          call exner_milieu(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
     223        endif
    220224        CALL massdair(p,masse)
    221225

  • trunk/libf/dyn3d/logic.h

    r119 r124  
    33!
    44!
    5 !
     5! NB: keep items of different kinds in seperate common blocs to avoid
     6!     "misaligned commons" issues
    67!-----------------------------------------------------------------------
    78! INCLUDE 'logic.h'
    89
    9       COMMON/logic/ purmats,iflag_phys,forward,leapf,apphys,            &
     10      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
    1011     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
    1112     &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
    12      &  ,ok_limit,ok_etat0,iflag_trac
     13     &  ,ok_limit,ok_etat0,hybrid
    1314
     15      COMMON/logici/ iflag_phys,iflag_trac
     16     
    1417      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
    1518     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
    1619     &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
    1720     &  ,ok_limit,ok_etat0
     21      logical hybrid ! vertcal coordinate is hybrid if true (sigma otherwise)
    1822
    1923      INTEGER iflag_phys,iflag_trac

  • trunk/libf/dyn3d/sortvarc.F

    r101 r124  
    6363
    6464c-----------------------------------------------------------------------
     65! Ehouarn: when no initialization fields from file, resetvarc should be
     66!          set to false
     67       if (firstcal) then
     68         if (.not.read_start) then
     69           resetvarc=.true.
     70         endif
     71       endif
    6572
    6673       dtvrs1j   = dtvr/daysec

  • trunk/libf/dyn3dpar/comvert.h

    r109 r124  
    77      COMMON/comvert/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),      &
    88     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
    9      &               aps(llm),bps(llm)
     9     &               aps(llm),bps(llm),scaleheight
    1010
    1111      REAL ap,bp,presnivs,dpres,pa,preff,nivsigs,nivsig,aps,bps
     12      real scaleheight ! atmospheric (reference) scale height (km)
    1213
    1314 !-----------------------------------------------------------------------

  • trunk/libf/dyn3dpar/disvert_noterre.F

    r110 r124  
    44c    Nouvelle version 100% Mars !!
    55c    On l'utilise aussi pour Venus et Titan, legerment modifiee.
     6
     7#ifdef CPP_IOIPSL
     8      use IOIPSL
     9#else
     10! if not using IOIPSL, we still need to use (a local version of) getin
     11      use ioipsl_getincom
     12#endif
    613
    714      IMPLICIT NONE
    … …  
    1219#include "comconst.h"
    1320#include "logic.h"
     21#include "iniprint.h"
    1422c
    1523c=======================================================================
    … …  
    2432      INTEGER l,ll
    2533      REAL snorm
    26       REAL alpha,beta,gama,delta,deltaz,h,quoi,quand
     34      REAL alpha,beta,gama,delta,deltaz
     35      real quoi,quand
    2736      REAL zsig(llm),sig(llm+1)
    2837      INTEGER np,ierr
    … …  
    4453c-----------------------------------------------------------------------
    4554c
     55! Initializations:
    4656      pi=2.*ASIN(1.)
     57     
     58      hybrid=.true. ! default value for hybrid (ie: use hybrid coordinates)
     59      CALL getin('hybrid',hybrid)
     60      write(lunout,*)'disvert_noterre: hybrid=',hybrid
    4761
    4862! Ouverture possible de fichiers typiquement E.T.
    … …  
    6781c        <-> energie cinetique, d'apres la note de Frederic Hourdin...
    6882
    69          PRINT*,'*****************************'
    70          PRINT*,'WARNING Lecture de esasig.def'
    71          PRINT*,'*****************************'
    72          READ(99,*) h
     83         write(lunout,*)'*****************************'
     84         write(lunout,*)'WARNING reading esasig.def'
     85         write(lunout,*)'*****************************'
     86         READ(99,*) scaleheight
    7387         READ(99,*) dz0
    7488         READ(99,*) dz1
    … …  
    7690         CLOSE(99)
    7791
    78          dz0=dz0/h
    79          dz1=dz1/h
     92         dz0=dz0/scaleheight
     93         dz1=dz1/scaleheight
    8094
    8195         sig1=(1.-dz1)/tanh(.5*(llm-1)/nhaut)
    … …  
    116130
    117131      ELSE IF(ierr4.eq.0) then
    118          PRINT*,'****************************'
    119          PRINT*,'Lecture de z2sig.def'
    120          PRINT*,'****************************'
    121 
    122          READ(99,*) h
     132         write(lunout,*)'****************************'
     133         write(lunout,*)'Reading z2sig.def'
     134         write(lunout,*)'****************************'
     135
     136         READ(99,*) scaleheight
    123137         do l=1,llm
    124138            read(99,*) zsig(l)
    … …  
    128142         sig(1) =1
    129143         do l=2,llm
    130            sig(l) = 0.5 * ( exp(-zsig(l)/h) + exp(-zsig(l-1)/h) )
     144           sig(l) = 0.5 * ( exp(-zsig(l)/scaleheight) +
     145     &                      exp(-zsig(l-1)/scaleheight) )
    131146         end do
    132147         sig(llm+1) =0
    … …  
    134149c-----------------------------------------------------------------------
    135150      ELSE
    136          write(*,*) 'didn t you forget something ??? '
    137          write(*,*) 'We need the file  z2sig.def ! (OR esasig.def) '
     151         write(lunout,*) 'didn t you forget something ??? '
     152         write(lunout,*) 'We need file  z2sig.def ! (OR esasig.def)'
    138153         stop
    139154      ENDIF
    … …  
    156171c-----------------------------------------------------------------------
    157172c
    158 c se trouvait dans Titan... Verifier...
    159 c     h = 40.
    160 
    161       if (1.eq.1) then  ! toujours hybrides
    162          write(*,*) "*******************************"
    163          write(*,*) "Systeme en coordonnees hybrides"
    164          write(*,*)
     173
     174      if (hybrid) then  ! use hybrid coordinates
     175         write(lunout,*) "*********************************"
     176         write(lunout,*) "Using hybrid vertical coordinates"
     177         write(lunout,*)
    165178c        Coordonnees hybrides avec mod
    166179         DO l = 1, llm
    … …  
    172185         bp(llmp1) = 0.
    173186         ap(llmp1) = 0.
    174       else
    175          write(*,*) "****************************"
    176          write(*,*) "Systeme en coordonnees sigma"
    177          write(*,*)
     187      else ! use sigma coordinates
     188         write(lunout,*) "********************************"
     189         write(lunout,*) "Using sigma vertical coordinates"
     190         write(lunout,*)
    178191c        Pour ne pas passer en coordonnees hybrides
    179192         DO l = 1, llm
    … …  
    186199      bp(llmp1) =   0.
    187200
    188       PRINT *,' BP '
    189       PRINT *,  bp
    190       PRINT *,' AP '
    191       PRINT *,  ap
     201      write(lunout,*)' BP '
     202      write(lunout,*)  bp
     203      write(lunout,*)' AP '
     204      write(lunout,*)  ap
    192205
    193206c     Calcul au milieu des couches :
    … …  
    197210c     (on met la meme distance (en log pression)  entre P(llm)
    198211c     et P(llm -1) qu'entre P(llm-1) et P(llm-2) ) est
    199 c     Specifique.  Ce choix est spécifié ici ET dans exner_hyb.F
     212c     Specifique.  Ce choix est spécifié ici ET dans exner_milieu.F
    200213
    201214      DO l = 1, llm-1
    … …  
    204217      ENDDO
    205218     
    206 c     if (hybrid) then
     219      if (hybrid) then
    207220         aps(llm) = aps(llm-1)**2 / aps(llm-2)
    208221         bps(llm) = 0.5*(bp(llm) + bp(llm+1))
    209 c     else
    210 c        bps(llm) = bps(llm-1)**2 / bps(llm-2)
    211 c        aps(llm) = 0.
    212 c     end if
    213 
    214       PRINT *,' BPs '
    215       PRINT *,  bps
    216       PRINT *,' APs'
    217       PRINT *,  aps
     222      else
     223         bps(llm) = bps(llm-1)**2 / bps(llm-2)
     224         aps(llm) = 0.
     225      end if
     226
     227      write(lunout,*)' BPs '
     228      write(lunout,*)  bps
     229      write(lunout,*)' APs'
     230      write(lunout,*)  aps
    218231
    219232      DO l = 1, llm
    220233       presnivs(l) = aps(l)+bps(l)*preff
    221        pseudoalt(l) = -h*log(presnivs(l)/preff)
     234       pseudoalt(l) = -scaleheight*log(presnivs(l)/preff)
    222235      ENDDO
    223236
    224       PRINT *,' PRESNIVS'
    225       PRINT *,presnivs
    226       PRINT *,'Pseudo altitude des Presnivs : (avec H = ',h,' km)'
    227       PRINT *,pseudoalt
     237      write(lunout,*)' PRESNIVS'
     238      write(lunout,*)presnivs
     239      write(lunout,*)'Pseudo altitude of Presnivs : (for a scale ',
     240     &                'height of ',scaleheight,' km)'
     241      write(lunout,*)pseudoalt
    228242
    229243c     --------------------------------------------------
    … …  
    232246c     --------------------------------------------------
    233247c     open (53,file='testhybrid.tab')
    234 c     h=15.5
     248c     scaleheight=15.5
    235249c     do iz=0,34
    236250c       z = -5 + min(iz,34-iz)
    237251c     approximation of scale height for Venus
    238 c       h = 15.5 - z/55.*10.
    239 c       ps = preff*exp(-z/h)
    240 c       zsig(1)= -h*log((aps(1) + bps(1)*ps)/preff)
     252c       scaleheight = 15.5 - z/55.*10.
     253c       ps = preff*exp(-z/scaleheight)
     254c       zsig(1)= -scaleheight*log((aps(1) + bps(1)*ps)/preff)
    241255c       do l=2,llm
    242256c     approximation of scale height for Venus
    243257c          if (zsig(l-1).le.55.) then
    244 c             h = 15.5 - zsig(l-1)/55.*10.
     258c             scaleheight = 15.5 - zsig(l-1)/55.*10.
    245259c          else
    246 c             h = 5.5 - (zsig(l-1)-55.)/35.*2.
     260c             scaleheight = 5.5 - (zsig(l-1)-55.)/35.*2.
    247261c          endif
    248 c          zsig(l)= zsig(l-1)-h*
     262c          zsig(l)= zsig(l-1)-scaleheight*
    249263c    .    log((aps(l) + bps(l)*ps)/(aps(l-1) + bps(l-1)*ps))
    250264c       end do

  • trunk/libf/dyn3dpar/exner_milieu.F

    r109 r124  
     1!
     2! $Id $
     3!
    14      SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
    25c
    … …  
    4447      REAL xpn, xps
    4548      REAL SSUM
    46       EXTERNAL filtreg, SSUM
     49      EXTERNAL SSUM
     50
     51      if (llm.eq.1) then
     52        ! Specific behaviour for Shallow Water (1 vertical layer) case
     53     
     54        ! Sanity checks
     55        if (kappa.ne.1) then
     56          call abort_gcm("exner_hyb",
     57     &    "kappa!=1 , but running in Shallow Water mode!!",42)
     58        endif
     59        if (cpp.ne.r) then
     60        call abort_gcm("exner_hyb",
     61     &    "cpp!=r , but running in Shallow Water mode!!",42)
     62        endif
     63       
     64        ! Compute pks(:),pk(:),pkf(:)
     65       
     66        DO   ij  = 1, ngrid
     67          pks(ij) = (cpp/preff) * ps(ij)
     68          pk(ij,1) = .5*pks(ij)
     69        ENDDO
     70       
     71        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
     72        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 )
     73       
     74        ! our work is done, exit routine
     75        return
     76      endif ! of if (llm.eq.1)
     77
    4778     
    4879c     -------------

  • trunk/libf/dyn3dpar/exner_milieu_p.F

    r109 r124  
    1       SUBROUTINE  exner_milieu ( ngrid, ps, p,beta, pks, pk, pkf )
     1!
     2! $Id $
     3!
     4      SUBROUTINE  exner_milieu_p ( ngrid, ps, p,beta, pks, pk, pkf )
    25c
    36c     Auteurs :  F. Forget , Y. Wanherdrick
    … …  
    5053c$OMP BARRIER
    5154
     55      if (llm.eq.1) then
     56        ! Specific behaviour for Shallow Water (1 vertical layer) case
     57     
     58        ! Sanity checks
     59        if (kappa.ne.1) then
     60          call abort_gcm("exner_hyb",
     61     &    "kappa!=1 , but running in Shallow Water mode!!",42)
     62        endif
     63        if (cpp.ne.r) then
     64        call abort_gcm("exner_hyb",
     65     &    "cpp!=r , but running in Shallow Water mode!!",42)
     66        endif
     67       
     68        ! Compute pks(:),pk(:),pkf(:)
     69        ijb=ij_begin
     70        ije=ij_end
     71!$OMP DO SCHEDULE(STATIC)
     72        DO ij=ijb, ije
     73          pks(ij)=(cpp/preff)*ps(ij)
     74          pk(ij,1) = .5*pks(ij)
     75          pkf(ij,1)=pk(ij,1)
     76        ENDDO
     77!$OMP ENDDO
     78
     79!$OMP MASTER
     80      if (pole_nord) then
     81        DO  ij   = 1, iim
     82          ppn(ij) = aire(   ij   ) * pks(  ij     )
     83        ENDDO
     84        xpn      = SSUM(iim,ppn,1) /apoln
     85 
     86        DO ij   = 1, iip1
     87          pks(   ij     )  =  xpn
     88          pk(ij,1) = .5*pks(ij)
     89          pkf(ij,1)=pk(ij,1)
     90        ENDDO
     91      endif
     92     
     93      if (pole_sud) then
     94        DO  ij   = 1, iim
     95          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
     96        ENDDO
     97        xps      = SSUM(iim,pps,1) /apols
     98 
     99        DO ij   = 1, iip1
     100          pks( ij+ip1jm )  =  xps
     101          pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
     102          pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
     103        ENDDO
     104      endif
     105!$OMP END MASTER
     106
     107        jjb=jj_begin
     108        jje=jj_end
     109        CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
     110
     111        ! our work is done, exit routine
     112        return
     113      endif ! of if (llm.eq.1)
     114
    52115c     -------------
    53116c     Calcul de pks

  • trunk/libf/dyn3dpar/gcm.F

    r101 r124  
    487487c   --------------------------------------
    488488      IF (config_inca /= 'none') THEN
     489#ifdef INCA
    489490!$OMP PARALLEL
    490 #ifdef INCA
    491491         CALL init_inca_dim(klon_omp,llm,iim,jjm,
    492492     $        rlonu,rlatu,rlonv,rlatv)
    493 #endif
    494493!$OMP END PARALLEL
     494#endif
    495495      END IF
    496496
    … …  
    574574c       write(78,*) 'q',q
    575575
    576 c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logic/)
     576c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
    577577      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,
    578578     .              time_0)

  • trunk/libf/dyn3dpar/iniacademic.F90

    r66 r124  
    207207        enddo
    208208     enddo
    209      write(lunout,*) 'Iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
     209!     write(lunout,*) 'Iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
    210210
    211211     ! 3. Initialize fields (if necessary)
    … …  
    217217
    218218        CALL pression ( ip1jmp1, ap, bp, ps, p       )
    219         CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
     219        if (planet_type=="earth") then
     220          CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
     221        else
     222          call exner_milieu(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
     223        endif
    220224        CALL massdair(p,masse)
    221225

  • trunk/libf/dyn3dpar/leapfrog_p.F

    r110 r124  
    600600!$OMP END DO
    601601
     602c$OMP MASTER
    602603      IF (prt_level>9) THEN
    603604        WRITE(lunout,*)"leapfrog_p: Iteration No",True_itau
    … …  
    15201521c$OMP MASTER
    15211522              nbetat = nbetatdem
     1523c$OMP END MASTER
     1524c$OMP BARRIER
    15221525
    15231526! ADAPTATION GCM POUR CP(T)
    … …  
    15301533      enddo
    15311534!$OMP END DO
     1535
     1536!$OMP MASTER
    15321537!              CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
    15331538      CALL geopot_p  ( ip1jmp1, tsurpk  , pk , pks,  phis  , phi   )
    … …  
    15991604                  endif
    16001605              endif ! of if (output_grads_dyn)
    1601 c$OMP END MASTER
     1606!$OMP END MASTER
    16021607             endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
    16031608            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
    … …  
    17441749c$OMP MASTER
    17451750                nbetat = nbetatdem
     1751c$OMP END MASTER
     1752c$OMP BARRIER
     1753
    17461754! ADAPTATION GCM POUR CP(T)
    17471755                call tpot2t_p(ip1jmp1,llm,teta,temp,pk)
    … …  
    17541762                enddo
    17551763!$OMP END DO
     1764
     1765!$OMP MASTER
    17561766!                CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
    17571767                CALL geopot_p(ip1jmp1,tsurpk,pk,pks,phis,phi)
    … …  
    18181828                endif ! of if (output_grads_dyn)
    18191829
    1820 c$OMP END MASTER
     1830!$OMP END MASTER
    18211831              ENDIF ! of IF(MOD(itau,iecri).EQ.0)
    18221832

  • trunk/libf/dyn3dpar/logic.h

    r119 r124  
    33!
    44!
    5 !
     5! NB: keep items of different kinds in seperate common blocs to avoid
     6!     "misaligned commons" issues
    67!-----------------------------------------------------------------------
    78! INCLUDE 'logic.h'
    89
    9       COMMON/logic/ purmats,iflag_phys,forward,leapf,apphys,            &
     10      COMMON/logicl/ purmats,forward,leapf,apphys,                      &
    1011     &  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus       &
    1112     &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
    12      &  ,ok_limit,ok_etat0,iflag_trac
     13     &  ,ok_limit,ok_etat0,hybrid
    1314
     15      COMMON/logici/ iflag_phys,iflag_trac
     16     
    1417      LOGICAL purmats,forward,leapf,apphys,statcl,conser,               &
    1518     & apdiss,apdelq,saison,ecripar,fxyhypb,ysinus                      &
    1619     &  ,read_start,ok_guide,ok_strato,ok_gradsfile                     &
    1720     &  ,ok_limit,ok_etat0
     21      logical hybrid ! vertcal coordinate is hybrid if true (sigma otherwise)
    1822
    1923      INTEGER iflag_phys,iflag_trac
    20 !$OMP THREADPRIVATE(/logic/)
     24!$OMP THREADPRIVATE(/logicl/)
     25!$OMP THREADPRIVATE(/logici/)
    2126!-----------------------------------------------------------------------

  • trunk/libf/dyn3dpar/sortvarc.F

    r101 r124  
    6363
    6464c-----------------------------------------------------------------------
     65! Ehouarn: when no initialization fields from file, resetvarc should be
     66!          set to false
     67       if (firstcal) then
     68         if (.not.read_start) then
     69           resetvarc=.true.
     70         endif
     71       endif
    6572
    6673       dtvrs1j   = dtvr/daysec
    … …  
    121128     *               cosphi(ij)
    122129          ENDDO
    123           angl(l) = radsg *
     130          angl(l) = rad *
    124131     s    (SSUM(ip1jm-iip1,ge(iip2),1)-SSUM(jjm-1,ge(iip2),iip1))
    125132      ENDDO

  • trunk/libf/dyn3dpar/top_bound_p.F

    r110 r124  
    5050!     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/)
    5151      LOGICAL,SAVE :: first=.true.
     52     
     53      REAL zkm
    5254      INTEGER j,l,jjb,jje
    5355
    … …  
    100102      IF (pole_sud) jje=jj_end-1
    101103
    102 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    103104      if (mode_top_bound.ge.2) then
     105!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    104106       do l=1,llm
    105107        do j=jjb,jje
    … …  
    116118        enddo
    117119       enddo
     120!$OMP END DO NOWAIT   
    118121      else
     122!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    119123       do l=1,llm
    120124        do j=jjb,jje
    … …  
    122126        enddo
    123127       enddo
    124       endif
    125 c$OMP END DO NOWAIT   
     128!$OMP END DO NOWAIT
     129      endif
    126130
    127131C   AMORTISSEMENTS LINEAIRES:
    128132
    129 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    130133      if (mode_top_bound.ge.1) then
     134!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    131135       do l=1,llm
    132136        do j=jjb,jje
    … …  
    136140        enddo
    137141       enddo
    138       endif
    139 c$OMP END DO NOWAIT
     142!$OMP END DO NOWAIT
     143      endif
    140144
    141145C POUR U ET H
    … …  
    148152      IF (pole_sud)  jje=jj_end-1
    149153
    150 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    151154      if (mode_top_bound.ge.2) then
     155!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    152156       do l=1,llm
    153157        do j=jjb,jje
    … …  
    161165        enddo
    162166       enddo
     167!$OMP END DO NOWAIT
    163168      else
     169!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    164170       do l=1,llm
    165171        do j=jjb,jje
    … …  
    167173        enddo
    168174       enddo
    169       endif
    170 c$OMP END DO NOWAIT
    171 
    172 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
     175!$OMP END DO NOWAIT
     176      endif
     177
    173178      if (mode_top_bound.ge.3) then
     179!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    174180       do l=1,llm
    175181        do j=jjb,jje
    … …  
    183189        enddo
    184190       enddo
    185       endif
    186 c$OMP END DO NOWAIT
     191!$OMP END DO NOWAIT
     192      endif
    187193
    188194C   AMORTISSEMENTS LINEAIRES:
    189195
    190 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    191196      if (mode_top_bound.ge.1) then
     197!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    192198       do l=1,llm
    193199        do j=jjb,jje
    … …  
    197203        enddo
    198204       enddo
    199       endif
    200 c$OMP END DO NOWAIT
     205!$OMP END DO NOWAIT
     206      endif
    201207     
    202 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    203208      if (mode_top_bound.ge.3) then
     209!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    204210       do l=1,llm
    205211        do j=jjb,jje
    … …  
    209215        enddo
    210216       enddo
    211       endif
    212 c$OMP END DO NOWAIT
    213      
    214 
     217!$OMP END DO NOWAIT
     218      endif     
     219
     220!$OMP BARRIER
    215221      RETURN
    216222      END
Note: See TracChangeset for help on using the changeset viewer.