Changeset 124 for trunk/libf/dyn3d


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/libf/dyn3d
Files:
6 edited

Legend:

Unmodified
Added
Removed

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