Changeset 124 for trunk/libf/dyn3dpar


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/dyn3dpar
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • 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.