Changeset 108 for trunk/libf/dyn3d


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

Sebastien Lebonnois: sponge layer et dissip horizontale.

Location:
trunk/libf/dyn3d
Files:
6 edited

Legend:

Unmodified
Added
Removed

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