Changeset 108


Ignore:
Timestamp:
Apr 12, 2011, 11:16:02 AM (14 years ago)
Author:
slebonnois
Message:

Sebastien Lebonnois: sponge layer et dissip horizontale.

Location:
trunk
Files:
4 added
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/chantiers/commit_importants.log

    r105 r108  
    736736Avec la v105, la compilation Titan se fait sans problemes.
    737737
    738 
     738**-----------
     739(commit 106: ajustements dans phytitan)
     740**-----------
     741- Sorties des opacites (*_hist*.h)
     742- g(z) pris en compte (physiq.F, zzlay)
     743- correction topo (physiq.F, zzlev(i,1))
     744- correction lell.F et lell_light.F (et aussi radtitan.F du coup)
     745- remise a jour des sorties
     746- suphec.F : RDAY=1.37889e6 et non RSIDAY
     747**-----------
     748
     749*********************
     750**** commit_v108 ****
     751*********************
     752
     753Modifs Sebastien Lebonnois dans dyn3d[par] pour sponge layer et dissip horizontale.
     754
     755** Correction bug
     756**---------------
     757
     758* Modif de dyn3d[par]/calfis[_p].F
     759  => ztfi doublement incremente...
     760
     761** Sponge layer
     762**-------------
     763
     764* Modif de dyn3d[par]/leapfrog[_p].F
     765  Utilisation de d*top dans appel a top_bound.F
     766  Initialisation des d* (dis, fi, top)
     767
     768* Modif de dyn3d[par]/top_bound[_p].F
     769  Implementation de mode_top_bound
     770  Pas d'incrementation
     771
     772* Modif de dyn3d[par]/conf_gcm.F
     773  Parametres de top_bound
     774
     775* Modif de dyn3d[par]/comconst.h
     776  Parametres de top_bound
     777
     778* Modif de deftank/[titan/venus]/gcm.def
     779  Ajout ok_strato
     780  Ajout parametres de top_bound
     781
     782* Ecriture doc
     783  Ajout top_bound.tex (et .pdf) dans la documentation.
     784
     785** Dissipation horizontale
     786**------------------------
     787
     788* Modif de dyn3d[par]/inidissip.F
     789  passage a dissip_fac_mid automatique
     790  passage a dissip_fac_up sur modele martien si ok_strato
     791
     792* Modif de dyn3d[par]/conf_gcm.F
     793  Parametres dissip_*
     794
     795* Modif de dyn3d[par]/comconst.h
     796  Parametres dissip_*
     797
     798* Modif de deftank/[titan/venus]/gcm.def
     799  Ajout parametres dissip_*   !!! REGLAGES A FAIRE !!!
     800
     801* Ecriture doc
     802  Ajout dissip_horiz.tex (et .pdf) dans la documentation.
     803
     804Reste la discretisation verticale a voir.
     805

  • trunk/deftanks/titan/gcm.def

    r6 r108  
    4141## avec ou sans traceurs                                                 
    4242iflag_trac=1
     43##  Avec ou sans strato // i.e. Couche eponge et second palier pour dissip horiz.
     44ok_strato=y
     45## Dissipation horizontale
     46dissip_fac_mid=2.
     47dissip_fac_up=10.
     48# deltaz et hdelta en km
     49dissip_deltaz=10.
     50dissip_hdelta=5.
     51# pupstart en Pa
     52dissip_pupstart=1.e3
     53## Couche eponge
     54#   1: dans les 4 derniers niveaux
     55#   2: dans les couches de pression plus faible que 100 fois la pression de la derniere couche
     56iflag_top_bound=1
     57## Mode Couche eponge
     58#   mode = 0 : pas de sponge
     59#   mode = 1 : u et v -> 0
     60#   mode = 2 : u et v -> moyenne zonale
     61#   mode = 3 : u, v et h -> moyenne zonale
     62mode_top_bound=1
     63#  Coefficient pour la couche eponge (valeur derniere couche)
     64tau_top_bound=4.e-5
    4365#
    4466## longitude en degres du centre du zoom                                 

  • trunk/deftanks/venus/gcm.def

    r6 r108  
    4141## avec ou sans traceurs                                                 
    4242iflag_trac=0
     43##  Avec ou sans strato // i.e. Couche eponge et second palier pour dissip horiz.
     44ok_strato=y
     45## Dissipation horizontale
     46dissip_fac_mid=2.
     47dissip_fac_up=2.
     48# deltaz et hdelta en km
     49dissip_deltaz=10.
     50dissip_hdelta=5.
     51# pupstart en Pa
     52dissip_pupstart=1.e3
     53## Couche eponge
     54#   1: dans les 4 derniers niveaux
     55#   2: dans les couches de pression plus faible que 100 fois la pression de la derniere couche
     56iflag_top_bound=1
     57## Mode Couche eponge
     58#   mode = 0 : pas de sponge
     59#   mode = 1 : u et v -> 0
     60#   mode = 2 : u et v -> moyenne zonale
     61#   mode = 3 : u, v et h -> moyenne zonale
     62mode_top_bound=1
     63#  Coefficient pour la couche eponge (valeur derniere couche)
     64tau_top_bound=1.e-5
    4365#
    4466## longitude en degres du centre du zoom                                 

  • trunk/documentation/cpdet.tex

    r5 r108  
    112112\section{Pratical aspects in the code}
    113113
    114 A specific file has been added to the dynamical core, cpdet.F, which includes
     114A specific file has been added to the dynamical core, \textsf{cpdet.F}, which includes
    115115all the needed routines to take the $C_p(T)$ possibility into account.
    116116These routines take advantage of the keyword \textsf{planet\_type} to

  • trunk/libf/dyn3d/calfis.F

    r97 r108  
    587587
    588588! ADAPTATION GCM POUR CP(T)
    589          ztfi=ztfi+zdtfi*dtphys
    590589      call t2tpot(ngridmx*llm,ztfi,zteta,zpk)
    591590

  • trunk/libf/dyn3d/comconst.h

    r53 r108  
    66
    77      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
    8      &                 iflag_top_bound
     8     &                 iflag_top_bound,mode_top_bound
    99      COMMON/comconstr/dtvr,daysec,                                     &
    1010     & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
    11      &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
    12      &                   ,tau_top_bound,                                &
     11     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
     12     & ,dissip_pupstart  ,tau_top_bound,                                &
    1313     & daylen,year_day,molmass, ihf
    1414      COMMON/cpdetvenus/nu_venus,t0_venus
    … …  
    2828      REAL g ! (m/s2) gravity
    2929      REAL omeg ! (rad/s) rotation rate of the planet
    30       REAL dissip_factz,dissip_deltaz,dissip_zref
    31       INTEGER iflag_top_bound
     30      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
     31      REAL dissip_pupstart
     32      INTEGER iflag_top_bound,mode_top_bound
    3233      REAL tau_top_bound
    3334      REAL daylen ! length of solar day, in 'standard' day length

  • trunk/libf/dyn3d/conf_gcm.F

    r97 r108  
    303303! Parametres controlant la variation sur la verticale des constantes de
    304304! dissipation.
    305 ! Pour le moment actifs uniquement dans la version a 39 niveaux
    306 ! avec ok_strato=y
    307 
    308        dissip_factz=4.
    309        dissip_deltaz=10.
    310        dissip_zref=30.
    311        CALL getin('dissip_factz',dissip_factz )
     305! Actifs uniquement avec ok_strato=y
     306
     307       dissip_fac_mid=2.
     308       dissip_fac_up=10.
     309       dissip_deltaz=10.! Intervalle (km) pour le changement mid / up
     310       dissip_hdelta=5. ! scale height (km) dans la zone de la transition(m)
     311       dissip_pupstart=1.e3  ! pression (Pa) au bas la transition mid / up
     312       CALL getin('dissip_fac_mid',dissip_fac_mid )
     313       CALL getin('dissip_fac_up',dissip_fac_up )
    312314       CALL getin('dissip_deltaz',dissip_deltaz )
    313        CALL getin('dissip_zref',dissip_zref )
    314 
     315       CALL getin('dissip_hdelta',dissip_hdelta )
     316       CALL getin('dissip_pupstart',dissip_pupstart )
     317
     318! Parametres controlant la couche eponge
     319! Actifs uniquement avec ok_strato=y
     320!   mode = 0 : pas de sponge
     321!   mode = 1 : u et v -> 0
     322!   mode = 2 : u et v -> moyenne zonale
     323!   mode = 3 : u, v et h -> moyenne zonale
    315324       iflag_top_bound=1
     325       mode_top_bound=3
    316326       tau_top_bound=1.e-5
    317327       CALL getin('iflag_top_bound',iflag_top_bound)
     328       CALL getin('mode_top_bound',mode_top_bound)
    318329       CALL getin('tau_top_bound',tau_top_bound)
    319330

  • trunk/libf/dyn3d/inidissip.F

    r1 r108  
    3131      INTEGER l,ij,idum,ii
    3232      REAL tetamin
    33       REAL pseudoz
     33      REAL Pup
    3434
    3535      REAL ran1
    … …  
    177177c   --------------------------------------------------
    178178
    179       if (ok_strato .and. llm==39) then
    180          do l=1,llm
    181             pseudoz=8.*log(preff/presnivs(l))
    182             zvert(l)=1+
    183      s      (tanh((pseudoz-dissip_zref)/dissip_deltaz)+1.)/2.
    184      s      *(dissip_factz-1.)
    185          enddo
    186       else
    187          DO l=1,llm
    188             zvert(l)=1.
    189          ENDDO
    190          fact=2.
    191          DO l = 1, llm
    192             zz      = 1. - preff/presnivs(l)
    193             zvert(l)= fact -( fact-1.)/( 1.+zz*zz )
    194          ENDDO
     179c First step: going from 1 to dissip_fac_mid (in gcm.def)
     180c============
     181      DO l=1,llm
     182         zz      = 1. - preff/presnivs(l)
     183         zvert(l)= dissip_fac_mid -( dissip_fac_mid-1.)/( 1.+zz*zz )
     184      ENDDO
     185
     186         write(*,*) 'Dissipation : '
     187         write(*,*) 'Multiplication de la dissipation en altitude :',
     188         write(*,*) '  dissip_fac_mid =', dissip_fac_mid
     189
     190c Second step if ok_strato:  from dissip_fac_mid to dissip_fac_up (in gcm.def)
     191c==========================
     192c Utilisation de la fonction tangente hyperbolique pour augmenter
     193c arbitrairement la dissipation et donc la stabilite du modele a
     194c partir d'une certaine altitude.
     195
     196c   Le facteur multiplicatif de basse atmosphere etant deja pris
     197c   en compte, il faut diviser le facteur multiplicatif de haute
     198c   atmosphere par celui-ci.
     199
     200      if (ok_strato) then
     201
     202       Pup = dissip_pupstart*exp(-0.5*dissip_deltaz/dissip_hdelta)
     203       do l=1,llm
     204         zvert(l)= zvert(l)*(1.0+( (dissip_fac_up/dissip_fac_mid-1)
     205     &   *(1-(0.5*(1+tanh(-6./dissip_deltaz*
     206     &    (-dissip_hdelta*log(presnivs(l)/Pup))  ))))  ))
     207       enddo
     208
     209         write(*,*) '  dissip_fac_up =', dissip_fac_up
     210         write(*,*) 'Transition mid /up:  Pupstart,delta =',
     211     &             dissip_pupstart,'Pa', dissip_deltaz , '(km)'
     212
    195213      endif
    196214

  • trunk/libf/dyn3d/leapfrog.F

    r101 r108  
    111111      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
    112112      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
    113 
    114       real :: duspg(ip1jmp1,llm) ! for bilan_dyn
    115113
    116114c   variables pour le fichier histoire
    … …  
    463461c      -------------------
    464462         IF (ok_strato) THEN
    465            CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
     463           CALL top_bound( vcov,ucov,teta,masse,dutop,dvtop,dtetatop)
     464c dqtop=0, dptop=0
     465           CALL addfi( dtphys, leapf, forward   ,
     466     $                  ucov, vcov, teta , q   ,ps ,
     467     $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
    466468         ENDIF
    467 c dqtop=0, dptop=0
    468        
    469 c      ajout des tendances physiques:
    470 c      ------------------------------
    471           CALL addfi( dtphys, leapf, forward   ,
    472      $                  ucov, vcov, teta , q   ,ps ,
    473      $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
    474 c
     469
    475470c  Diagnostique de conservation de l'énergie : difference
    476471         IF (ip_ebil_dyn.ge.1 ) THEN
    … …  
    512507        ! Sponge layer (if any)
    513508        IF (ok_strato) THEN
    514           dutop(:,:)=0.
    515           dvtop(:,:)=0.
    516           dtetatop(:,:)=0.
    517           dqtop(:,:,:)=0.
    518           dptop(:)=0.
    519509          CALL top_bound(vcov,ucov,teta,masse,dutop,dvtop,dtetatop)
     510c dqtop=0, dptop=0
    520511          CALL addfi( dtvr, leapf, forward   ,
    521512     $                  ucov, vcov, teta , q   ,ps ,
    … …  
    669660                 CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
    670661     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
    671      &                 du,dudis,duspg,dufi)
     662     &                 du,dudis,dutop,dufi)
    672663#endif
    673664               END IF
    … …  
    801792                 CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
    802793     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
    803      &                 du,dudis,duspg,dufi)
     794     &                 du,dudis,dutop,dufi)
    804795#endif
    805796               ENDIF

  • trunk/libf/dyn3d/top_bound.F

    r1 r108  
    5252      LOGICAL,SAVE :: first=.true.
    5353
     54      REAL zkm
    5455      INTEGER j,l
    5556
    … …  
    6869             rdamp(llm-3)=tau_top_bound/8.
    6970         else if (iflag_top_bound.eq.2) then
    70 ! couce eponge dans toutes les couches de pression plus faible que
     71! couche eponge dans toutes les couches de pression plus faible que
    7172! 100 fois la pression de la derniere couche
    7273             rdamp(:)=tau_top_bound
    … …  
    7475         endif
    7576         first=.false.
    76          print*,'TOP_BOUND rdamp=',rdamp
     77         print*,'TOP_BOUND mode',mode_top_bound
     78         print*,'Coeffs pour la couche eponge a l equateur'
     79         print*,'p (Pa)  z(km)  tau (s)   dt*rdamp'
     80         do l=1,llm
     81           if (rdamp(l).ne.0.) then
     82            zkm        = phi(iip1/2,jjp1/2,l)/(1000*g)
     83          print*,presnivs(l),zkm,
     84     .          1./rdamp(l),
     85     .          dt*rdamp(l)
     86           endif
     87         enddo
    7788      endif
    7889
    7990      CALL massbar(masse,massebx,masseby)
    8091
    81       do l=1,llm
     92c   mode = 0 : pas de sponge
     93c   mode = 1 : u et v -> 0
     94c   mode = 2 : u et v -> moyenne zonale
     95c   mode = 3 : u, v et h -> moyenne zonale
     96
     97      if (mode_top_bound.ge.2) then
     98       do l=1,llm
    8299        do j=1,jjm
    83100          vzon(j,l)=0.
    … …  
    92109          vzon(j,l)=vzon(j,l)/zm
    93110        enddo
    94       enddo
     111       enddo
    95112
    96       do l=1,llm
    97         do i=1,iip1
    98           do j=1,jjm
    99             dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
    100           enddo
    101         enddo
    102       enddo
    103 
    104       do l=1,llm
     113       do l=1,llm
    105114        do j=2,jjm
    106115          uzon(j,l)=0.
    … …  
    112121          uzon(j,l)=uzon(j,l)/zm
    113122        enddo
    114       enddo
     123       enddo
     124      else
     125          vzon(:,:)=0.
     126          uzon(:,:)=0.
     127      endif
    115128
    116       do l=1,llm
     129      if (mode_top_bound.ge.3) then
     130       do l=1,llm
    117131        do j=2,jjm
    118132          zm=0.
    … …  
    124138          tzon(j,l)=tzon(j,l)/zm
    125139        enddo
    126       enddo
     140       enddo
     141      endif
    127142
    128143C   AMORTISSEMENTS LINEAIRES:
    129144
    130       do l=1,llm
     145      if (mode_top_bound.ge.1) then
     146       do l=1,llm
     147        do i=1,iip1
     148          do j=1,jjm
     149            dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
     150          enddo
     151        enddo
     152       enddo
     153
     154       do l=1,llm
    131155        do i=1,iip1
    132156          do j=2,jjm
    133             du(i,j,l)=du(i,j,l)
    134      s               -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
    135             dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
     157            du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
    136158          enddo
    137159        enddo
    138       enddo
     160       enddo
     161      endif
     162
     163      if (mode_top_bound.ge.3) then
     164       do l=1,llm
     165        do i=1,iip1
     166          do j=2,jjm
     167            dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l))
     168          enddo
     169        enddo
     170       enddo
     171      endif
    139172     
    140 
    141173      RETURN
    142174      END

  • trunk/libf/dyn3dpar/addfi_p.F

    r7 r108  
    2828c      pdufi(ip1jmp1,llm)     |
    2929c      pdvfi(ip1jm,llm)       |   respective
    30 c      pdhfi(ip1jmp1)         |      tendencies
     30c      pdhfi(ip1jmp1)         |      tendencies  (unit/s)
    3131c      pdtsfi(ip1jmp1)        |
    3232c

  • trunk/libf/dyn3dpar/calfis_p.F

    r97 r108  
    952952     
    953953! ADAPTATION GCM POUR CP(T)
    954 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    955       DO l=1,llm
    956         ztfi(1:klon,l)=ztfi(1:klon,l)+zdtfi(1:klon,l)*dtphys
    957       ENDDO
    958 !$OMP END DO
    959954      call t2tpot_p(ngridmx,llm,ztfi,zteta,zpk)
    960955

  • trunk/libf/dyn3dpar/comconst.h

    r52 r108  
    66
    77      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
    8      &                 iflag_top_bound
     8     &                 iflag_top_bound,mode_top_bound
    99      COMMON/comconstr/dtvr,daysec,                                     &
    1010     & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
    11      &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
    12      &                   ,tau_top_bound,                                &
     11     & ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta        &
     12     & ,dissip_pupstart  ,tau_top_bound,                                &
    1313     & daylen,year_day,molmass, ihf
    1414      COMMON/cpdetvenus/nu_venus,t0_venus
    … …  
    2828      REAL g ! (m/s2) gravity
    2929      REAL omeg ! (rad/s) rotation rate of the planet
    30       REAL dissip_factz,dissip_deltaz,dissip_zref
    31       INTEGER iflag_top_bound
     30      REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
     31      REAL dissip_pupstart
     32      INTEGER iflag_top_bound,mode_top_bound
    3233      REAL tau_top_bound
    3334      REAL daylen ! length of solar day, in 'standard' day length

  • trunk/libf/dyn3dpar/conf_gcm.F

    r97 r108  
    314314! Parametres controlant la variation sur la verticale des constantes de
    315315! dissipation.
    316 ! Pour le moment actifs uniquement dans la version a 39 niveaux
    317 ! avec ok_strato=y
    318 
    319        dissip_factz=4.
    320        dissip_deltaz=10.
    321        dissip_zref=30.
    322        CALL getin('dissip_factz',dissip_factz )
     316! Actifs uniquement avec ok_strato=y
     317
     318       dissip_fac_mid=2.
     319       dissip_fac_up=10.
     320       dissip_deltaz=10.! Intervalle (km) pour le changement mid / up
     321       dissip_hdelta=5. ! scale height (km) dans la zone de la transition(m)
     322       dissip_pupstart=1.e3  ! pression (Pa) au bas la transition mid / up
     323       CALL getin('dissip_fac_mid',dissip_fac_mid )
     324       CALL getin('dissip_fac_up',dissip_fac_up )
    323325       CALL getin('dissip_deltaz',dissip_deltaz )
    324        CALL getin('dissip_zref',dissip_zref )
    325 
     326       CALL getin('dissip_hdelta',dissip_hdelta )
     327       CALL getin('dissip_pupstart',dissip_pupstart )
     328
     329! Parametres controlant la couche eponge
     330! Actifs uniquement avec ok_strato=y
     331!   mode = 0 : pas de sponge
     332!   mode = 1 : u et v -> 0
     333!   mode = 2 : u et v -> moyenne zonale
     334!   mode = 3 : u, v et h -> moyenne zonale
    326335       iflag_top_bound=1
     336       mode_top_bound=3
    327337       tau_top_bound=1.e-5
    328338       CALL getin('iflag_top_bound',iflag_top_bound)
     339       CALL getin('mode_top_bound',mode_top_bound)
    329340       CALL getin('tau_top_bound',tau_top_bound)
    330341

  • trunk/libf/dyn3dpar/leapfrog_p.F

    r101 r108  
    113113      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
    114114
    115        !! Aymeric -- cp(T) comme dans leapfrog.F, SAVE OK ???
    116       REAL,SAVE :: duspg(ip1jmp1,llm) ! for bilan_dyn
    117 
     115c   tendances top_bound (sponge layer)
     116      REAL,SAVE :: dvtop(ip1jm,llm),dutop(ip1jmp1,llm)
     117      REAL,SAVE :: dtetatop(ip1jmp1,llm)
     118      REAL,SAVE :: dptop(ip1jmp1)
     119      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqtop
    118120
    119121c   variables pour le fichier histoire
    … …  
    242244            ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
    243245            ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
     246            ALLOCATE(dqtop(ip1jmp1,llm,nqtot))
    244247         END IF
    245248c$OMP END MASTER     
    … …  
    251254
    252255c$OMP MASTER
     256c INITIALISATIONS
     257        dudis(:,:)   =0.
     258        dvdis(:,:)   =0.
     259        dtetadis(:,:)=0.
     260        dutop(:,:)   =0.
     261        dvtop(:,:)   =0.
     262        dtetatop(:,:)=0.
     263        dqtop(:,:,:) =0.
     264        dptop(:)     =0.
     265        dufi(:,:)   =0.
     266        dvfi(:,:)   =0.
     267        dtetafi(:,:)=0.
     268        dqfi(:,:,:) =0.
     269        dpfi(:)     =0.
    253270      dq(:,:,:)=0.
     271
    254272      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    255273      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
    … …  
    914932c      ajout des tendances physiques:
    915933c      ------------------------------
    916          IF (ok_strato) THEN
    917            CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
    918          ENDIF
    919        
    920934          CALL addfi_p( dtphys, leapf, forward   ,
    921935     $                  ucov, vcov, teta , q   ,ps ,
    922936     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
    923937
     938c      Couche superieure :
     939c      -------------------
     940         IF (ok_strato) THEN
     941           CALL top_bound_p( vcov,ucov,teta,masse,dutop,dvtop,dtetatop)
     942           CALL addfi_p( dtphys, leapf, forward   ,
     943     $                  ucov, vcov, teta , q   ,ps ,
     944     $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
     945
     946         ENDIF
     947       
    924948c$OMP BARRIER
    925949c$OMP MASTER
    … …  
    14441468                 CALL bilan_dyn_p(dtvr*iperiod,dtvr*day_step*periodav,
    14451469     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
    1446      &                 du,dudis,duspg,dufi)
     1470     &                 du,dudis,dutop,dufi)
    14471471c$OMP END MASTER
    14481472              ENDIF !ok_dynzon
    … …  
    16701694                 CALL bilan_dyn_p(dtvr*iperiod,dtvr*day_step*periodav,
    16711695     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
    1672      &                 du,dudis,duspg,dufi)
     1696     &                 du,dudis,dutop,dufi)
    16731697
    16741698c$OMP END MASTER

  • trunk/libf/dyn3dpar/top_bound_p.F

    r1 r108  
    6464             rdamp(llm-3)=tau_top_bound/8.
    6565         else if (iflag_top_bound == 2) then
    66 ! couce eponge dans toutes les couches de pression plus faible que
     66! couche eponge dans toutes les couches de pression plus faible que
    6767! 100 fois la pression de la derniere couche
    6868             rdamp(:)=tau_top_bound
    … …  
    7070         endif
    7171         first=.false.
    72          print*,'TOP_BOUND rdamp=',rdamp
     72         print*,'TOP_BOUND mode',mode_top_bound
     73         print*,'Coeffs pour la couche eponge a l equateur'
     74         print*,'p (Pa)  z(km)  tau (s)   dt*rdamp'
     75         do l=1,llm
     76           if (rdamp(l).ne.0.) then
     77            zkm        = phi(iip1/2,jjp1/2,l)/(1000*g)
     78          print*,presnivs(l),zkm,
     79     .          1./rdamp(l),
     80     .          dt*rdamp(l)
     81           endif
     82         enddo
    7383c$OMP END MASTER
    7484c$OMP BARRIER
    … …  
    7787
    7888      CALL massbar_p(masse,massebx,masseby)
     89
     90c   mode = 0 : pas de sponge
     91c   mode = 1 : u et v -> 0
     92c   mode = 2 : u et v -> moyenne zonale
     93c   mode = 3 : u, v et h -> moyenne zonale
     94
     95C POUR V
     96
    7997C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
    8098
    … …  
    84102
    85103c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    86       do l=1,llm
     104      if (mode_top_bound.ge.2) then
     105       do l=1,llm
    87106        do j=jjb,jje
    88107          zm=0.
    … …  
    97116          vzon(j,l)=vzon(j,l)/zm
    98117        enddo
    99       enddo
     118       enddo
     119      else
     120       do l=1,llm
     121        do j=jjb,jje
     122          vzon(j,l)=0.
     123        enddo
     124       enddo
     125      endif
    100126c$OMP END DO NOWAIT   
    101127
    102 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    103       do l=1,llm
     128C   AMORTISSEMENTS LINEAIRES:
     129
     130c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     131      if (mode_top_bound.ge.1) then
     132       do l=1,llm
    104133        do j=jjb,jje
    105134          do i=1,iip1
    106             dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
    107           enddo
    108         enddo
    109       enddo
    110 c$OMP END DO NOWAIT
     135            dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
     136          enddo
     137        enddo
     138       enddo
     139      endif
     140c$OMP END DO NOWAIT
     141
     142C POUR U ET H
     143
     144C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
    111145
    112146      jjb=jj_begin
    … …  
    116150
    117151c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    118       do l=1,llm
     152      if (mode_top_bound.ge.2) then
     153       do l=1,llm
    119154        do j=jjb,jje
    120155          uzon(j,l)=0.
    … …  
    126161          uzon(j,l)=uzon(j,l)/zm
    127162        enddo
    128       enddo
     163       enddo
     164      else
     165       do l=1,llm
     166        do j=jjb,jje
     167          uzon(j,l)=0.
     168        enddo
     169       enddo
     170      endif
    129171c$OMP END DO NOWAIT
    130172
    131173c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
    132       do l=1,llm
     174      if (mode_top_bound.ge.3) then
     175       do l=1,llm
    133176        do j=jjb,jje
    134177          zm=0.
    … …  
    140183          tzon(j,l)=tzon(j,l)/zm
    141184        enddo
    142       enddo
     185       enddo
     186      endif
    143187c$OMP END DO NOWAIT
    144188
    … …  
    146190
    147191c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    148       do l=1,llm
     192      if (mode_top_bound.ge.1) then
     193       do l=1,llm
    149194        do j=jjb,jje
    150195          do i=1,iip1
    151             du(i,j,l)=du(i,j,l)
    152      s               -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
    153             dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
    154           enddo
    155        enddo
    156       enddo
     196            du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
     197          enddo
     198        enddo
     199       enddo
     200      endif
     201c$OMP END DO NOWAIT
     202     
     203c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     204      if (mode_top_bound.ge.3) then
     205       do l=1,llm
     206        do j=jjb,jje
     207          do i=1,iip1
     208            dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l))
     209          enddo
     210        enddo
     211       enddo
     212      endif
    157213c$OMP END DO NOWAIT
    158214     
Note: See TracChangeset for help on using the changeset viewer.