Ignore:
Timestamp:
Oct 11, 2007, 3:43:42 PM (17 years ago)
Author:
Laurent Fairhead
Message:

Mise a jour de la physique avec thermiques avec la version de FH d'aout 2007
LF

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/phytherm/thermcell_flux.F90

    r814 r852  
    2828
    2929      integer ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha
    30       integer ncorecfm4,ncorecfm5,ncorecfm6
     30      integer ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8
    3131     
    3232
    3333      REAL entr(ngrid,klev),detr(ngrid,klev)
    3434      REAL fm(ngrid,klev+1)
     35      REAL zfm
    3536
    3637      integer igout
    … …  
    4041      REAL f_old,ddd0,eee0,ddd,eee,zzz
    4142
    42       REAL fracmax
    43       save fracmax
    44 
    45       fracmax=0.5
     43      REAL fomass_max,alphamax
     44      save fomass_max,alphamax
     45
     46      fomass_max=0.5
     47      alphamax=0.7
    4648
    4749      ncorecfm1=0
    … …  
    5153      ncorecfm5=0
    5254      ncorecfm6=0
     55      ncorecfm7=0
     56      ncorecfm8=0
    5357      ncorecalpha=0
    5458
    … …  
    5660      fm(:,:)=0.
    5761     
     62      if (lev_out.ge.10) then
     63         write(lunout,*) 'Dans thermcell_flux 0'
     64         write(lunout,*) 'flux base ',f(igout)
     65         write(lunout,*) 'lmax ',lmax(igout)
     66         write(lunout,*) 'lalim ',lalim(igout)
     67         write(lunout,*) 'ig= ',igout
     68         write(lunout,*) ' l E*    A*     D*  '
     69         write(lunout,'(i4,3e15.5)') (l,entr_star(igout,l),alim_star(igout,l),detr_star(igout,l) &
     70     &    ,l=1,lmax(igout))
     71      endif
     72
    5873
    5974!-------------------------------------------------------------------------
    … …  
    6681            if (l.le.lmax(ig)) then
    6782               if (entr_star(ig,l).gt.1.) then
    68                     print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
     83                    print*,'WARNING thermcell_flux 1 ig,l,lmax(ig)',ig,l,lmax(ig)
     84                    print*,'entr_star(ig,l)',entr_star(ig,l)
     85                    print*,'alim_star(ig,l)',alim_star(ig,l)
     86                    print*,'detr_star(ig,l)',detr_star(ig,l)
     87!                   stop
     88               endif
     89            else
     90               if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0.) then
     91                    print*,'cas 1 : ig,l,lmax(ig)',ig,l,lmax(ig)
    6992                    print*,'entr_star(ig,l)',entr_star(ig,l)
    7093                    print*,'alim_star(ig,l)',alim_star(ig,l)
    … …  
    7295                    stop
    7396               endif
    74             else
    75                if (abs(entr_star(ig,l))+abs(alim_star(ig,l))+abs(detr_star(ig,l)).gt.0) then
    76                     print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
    77                     print*,'entr_star(ig,l)',entr_star(ig,l)
    78                     print*,'alim_star(ig,l)',alim_star(ig,l)
    79                     print*,'detr_star(ig,l)',detr_star(ig,l)
    80                     stop
    81                endif
    8297            endif
    8398         enddo
    … …  
    92107         detr(:,l)=f(:)*detr_star(:,l)
    93108      enddo
     109
     110      if (lev_out.ge.10) then
     111         write(lunout,*) 'Dans thermcell_flux 1'
     112         write(lunout,*) 'flux base ',f(igout)
     113         write(lunout,*) 'lmax ',lmax(igout)
     114         write(lunout,*) 'lalim ',lalim(igout)
     115         write(lunout,*) 'ig= ',igout
     116         write(lunout,*) ' l   E    D     W2'
     117         write(lunout,'(i4,3e15.5)') (l,entr(igout,l),detr(igout,l) &
     118     &    ,zw2(igout,l+1),l=1,lmax(igout))
     119      endif
    94120
    95121      fm(:,1)=0.
    … …  
    107133      enddo
    108134
    109       if (lev_out.ge.10) then
    110          write(lunout,*) 'Dans thermcell_flux 1'
    111          write(lunout,*) 'flux base ',f(igout)
    112          write(lunout,*) 'lmax ',lmax(igout)
    113          write(lunout,*) 'lalim ',lalim(igout)
    114          write(lunout,'(i6,i4,3e15.5)') (ig,l,entr(igout,l),detr(igout,l) &
    115      &    ,fm(igout,l+1),l=1,lmax(igout))
    116       endif
     135
     136
     137! Test provisoire : pour comprendre pourquoi on corrige plein de fois
     138! le cas fm6, on commence par regarder une premiere fois avant les
     139! autres corrections.
     140
     141      do l=1,klev
     142         do ig=1,ngrid
     143            if (detr(ig,l).gt.fm(ig,l)) then
     144               ncorecfm8=ncorecfm8+1
     145!              igout=ig
     146            endif
     147         enddo
     148      enddo
     149
     150      if (lev_out.ge.10) &
     151    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     152    &    ptimestep,masse,entr,detr,fm,'2  ')
    117153
    118154!-------------------------------------------------------------------------
    … …  
    130166         enddo
    131167      enddo
     168
     169      if (lev_out.ge.10) &
     170    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     171    &    ptimestep,masse,entr,detr,fm,'3  ')
    132172
    133173!-------------------------------------------------------------------------
    … …  
    151191      enddo
    152192
     193      if (lev_out.ge.10) &
     194    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     195    &    ptimestep,masse,entr,detr,fm,'4  ')
     196
    153197
    154198!-------------------------------------------------------------------------
    … …  
    167211      enddo
    168212
     213      if (lev_out.ge.10) &
     214    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     215    &    ptimestep,masse,entr,detr,fm,'5  ')
     216
    169217!-------------------------------------------------------------------------
    170218!detr ne peut pas etre superieur a fm
    … …  
    172220
    173221      if(1.eq.1) then
     222
    174223      do l=1,klev
    175224         do ig=1,ngrid
    … …  
    181230               ncorecfm6=ncorecfm6+1
    182231               detr(ig,l)=fm(ig,l)
    183                entr(ig,l)=fm(ig,l+1)
    184             endif
     232!              entr(ig,l)=fm(ig,l+1)
     233
     234! Dans le cas ou on est au dessus de la couche d'alimentation et que le
     235! detrainement est plus fort que le flux de masse, on stope le thermique.
     236               if (l.gt.lalim(ig)) then
     237                  lmax(ig)=l
     238                  fm(ig,l+1)=0.
     239                  entr(ig,l)=0.
     240               else
     241                  ncorecfm7=ncorecfm7+1
     242               endif
     243            endif
     244
     245            if(l.gt.lmax(ig)) then
     246               detr(ig,l)=0.
     247               fm(ig,l+1)=0.
     248               entr(ig,l)=0.
     249            endif
     250
    185251            if (entr(ig,l).lt.0.) then
    186252               print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
    … …  
    194260
    195261
     262      if (lev_out.ge.10) &
     263    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     264    &    ptimestep,masse,entr,detr,fm,'6  ')
     265
    196266!-------------------------------------------------------------------------
    197267!fm ne peut pas etre negatif
    … …  
    207277            endif
    208278            if (detr(ig,l).lt.0.) then
    209                print*,'ig,l,lmax(ig)',ig,l,lmax(ig)
     279               print*,'cas 2 : ig,l,lmax(ig)',ig,l,lmax(ig)
    210280               print*,'detr(ig,l)',detr(ig,l)
    211281               print*,'fm(ig,l)',fm(ig,l)
    … …  
    215285     enddo
    216286
     287      if (lev_out.ge.10) &
     288    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     289    &    ptimestep,masse,entr,detr,fm,'7  ')
     290
    217291!-----------------------------------------------------------------------
    218292!la fraction couverte ne peut pas etre superieure a 1           
    219293!-----------------------------------------------------------------------
     294
     295
     296!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     297! FH Partie a revisiter.
     298! Il semble qu'etaient codees ici deux optiques dans le cas
     299! F/ (rho *w) > 1
     300! soit limiter la hauteur du thermique en considerant que c'est
     301! la derniere chouche, soit limiter F a rho w.
     302! Dans le second cas, il faut en fait limiter a un peu moins
     303! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
     304! dans thermcell_main et qu'il semble de toutes facons deraisonable
     305! d'avoir des fractions de 1..
     306! Ci dessous, et dans l'etat actuel, le premier des  deux if est
     307! sans doute inutile.
     308!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     309
     310     if(1.eq.0) then
     311
    220312     do l=1,klev
    221313        do ig=1,ngrid
    222314           if (zw2(ig,l+1).gt.1.e-10) then
    223            if ((((fm(ig,l+1))/(rhobarz(ig,l+1)*zw2(ig,l+1))).gt.  &
    224      &      1.)) then
     315           zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax
     316           if ( fm(ig,l+1) .gt. zfm) then
    225317              f_old=fm(ig,l+1)
    226               fm(ig,l+1)=rhobarz(ig,l+1)*zw2(ig,l+1)
    227               zw2(ig,l+1)=0.
    228               zqla(ig,l+1)=0.
     318              fm(ig,l+1)=zfm
     319!             zw2(ig,l+1)=0.
     320!             zqla(ig,l+1)=0.
    229321              detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
    230               lmax(ig)=l+1
    231               zmax(ig)=zlev(ig,lmax(ig))
     322!             lmax(ig)=l+1
     323!             zmax(ig)=zlev(ig,lmax(ig))
    232324!             print*,'alpha>1',l+1,lmax(ig)
    233325              ncorecalpha=ncorecalpha+1
    … …  
    237329     enddo
    238330!
     331      if (lev_out.ge.10) &
     332    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     333    &    ptimestep,masse,entr,detr,fm,'8  ')
     334
     335     endif
    239336
    240337!-----------------------------------------------------------------------
    … …  
    248345            eee0=entr(ig,l)
    249346            ddd0=detr(ig,l)
    250             eee=entr(ig,l)-masse(ig,l)*fracmax/ptimestep
     347            eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep
    251348            ddd=detr(ig,l)-eee
    252349            if (eee.gt.0.) then
    … …  
    272369                         print*,'detr',detr(ig,l)
    273370                         print*,'masse',masse(ig,l)
    274                          print*,'fracmax',fracmax
    275                          print*,'masse(ig,l)*fracmax/ptimestep',masse(ig,l)*fracmax/ptimestep
     371                         print*,'fomass_max',fomass_max
     372                         print*,'masse(ig,l)*fomass_max/ptimestep',masse(ig,l)*fomass_max/ptimestep
    276373                         print*,'ptimestep',ptimestep
    277374                         print*,'lmax(ig)',lmax(ig)
    … …  
    309406    &     ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', &
    310407    &     ncorecfm6,'x fm6', &
     408    &     ncorecfm7,'x fm7', &
     409    &     ncorecfm8,'x fm8', &
    311410    &     ncorecalpha,'x alpha'
    312411      endif
    313412
     413      if (lev_out.ge.10) &
     414    &    call printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     415    &    ptimestep,masse,entr,detr,fm,'fin')
     416
    314417
    315418      return
    316419      end
     420
     421      subroutine printflux(ngrid,klev,lunout,igout,f,lmax,lalim, &
     422    &    ptimestep,masse,entr,detr,fm,descr)
     423
     424     implicit none
     425
     426      integer ngrid,klev,lunout,igout,l,lm
     427
     428      integer lmax(klev),lalim(klev)
     429      real ptimestep,masse(ngrid,klev),entr(ngrid,klev),detr(ngrid,klev)
     430      real fm(ngrid,klev+1),f(ngrid)
     431
     432      character*3 descr
     433
     434      lm=lmax(igout)+5
     435      if(lm.gt.klev) lm=klev
     436
     437      print*,'Impression jusque lm=',lm
     438
     439         write(lunout,*) 'Dans thermcell_flux '//descr
     440         write(lunout,*) 'flux base ',f(igout)
     441         write(lunout,*) 'lmax ',lmax(igout)
     442         write(lunout,*) 'lalim ',lalim(igout)
     443         write(lunout,*) 'ig= ',igout
     444         write(lunout,'(a3,4a14)') 'l','M','E','D','F'
     445         write(lunout,'(i4,4e14.4)') (l,masse(igout,l)/ptimestep, &
     446     &     entr(igout,l),detr(igout,l) &
     447     &    ,fm(igout,l+1),l=1,lm)
     448
     449
     450      do l=lmax(igout)+1,klev
     451          if (abs(entr(igout,l))+abs(detr(igout,l))+abs(fm(igout,l)).gt.0.) then
     452          print*,'cas 1 : igout,l,lmax(igout)',igout,l,lmax(igout)
     453          print*,'entr(igout,l)',entr(igout,l)
     454          print*,'detr(igout,l)',detr(igout,l)
     455          print*,'fm(igout,l)',fm(igout,l)
     456          stop
     457          endif
     458      enddo
     459
     460      return
     461      end
     462
Note: See TracChangeset for help on using the changeset viewer.