Changeset 4377


Ignore:
Timestamp:
Jan 8, 2023, 1:44:53 PM (2 years ago)
Author:
fhourdin
Message:

Modification replay et thermcell_down (atelier)

Location:
LMDZ6/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/thermcell_down.F90

    r4365 r4377  
    1    SUBROUTINE thermcell_updown_dq(ngrid,nlay,lmax,eup,dup,edn,ddn,masse,theta,dtheta,thetau,thetad)
    2 
    3 !--------------------------------------------------------------
    4 ! thermcell_updown_dq: calcul du transport d'un traceur en présence
     1   SUBROUTINE thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,eup,dup,edn,ddn,masse,theta)
     2
     3!--------------------------------------------------------------
     4! thermcell_updown_dq: calcul du transport d'un traceur en pr??sence
    55! d'up/down drafts
    66!--------------------------------------------------------------
    … …  
    1616
    1717   integer,intent(in) :: ngrid,nlay
     18   real,intent(in) :: ptimestep
    1819   real,intent(in), dimension(ngrid,nlay) :: eup,dup,edn,ddn,masse
    19    real,intent(in), dimension(ngrid,nlay) :: theta
    20    real,intent(out), dimension(ngrid,nlay) :: thetau,thetad,dtheta
    21    integer, intent(out), dimension(ngrid) :: lmax
     20   real,intent(inout), dimension(ngrid,nlay) :: theta
     21   integer, intent(in), dimension(ngrid) :: lmax
    2222
    2323   
    2424! Local
    2525
    26    real, dimension(ngrid,nlay+1) :: fup,fdn
     26   real, dimension(ngrid,nlay+1) :: fup,fdn,fthu,fthd,fthe,fthtot
     27   real, dimension(ngrid,nlay) :: thetau,thetad,dtheta
     28   real :: www
    2729
    2830   integer ig,ilay
    2931
    3032   fdn(:,:)=0.
     33   fup(:,:)=0.
     34   fthu(:,:)=0.
     35   fthd(:,:)=0.
     36   fthe(:,:)=0.
     37   fthtot(:,:)=0.
    3138   thetad(:,:)=0.
     39   thetau(:,:)=0.
    3240
    3341   ! lmax : indice tel que fu(kmax+1)=0
    … …  
    3543   ! Dans ce cas, pas besoin d'initialiser thetad(lmax) ( =theta(lmax) )
    3644
     45   print*,'ON PASSE BIEN PAR LA CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC'
     46   ! Boucle pour le downdraft
    3747   do ilay=nlay,1,-1
    3848      do ig=1,ngrid
    39          if (ilay.le.lmax(ig)) then
     49         if (ilay.le.lmax(ig) .and. lmax(ig)>1) then
    4050            fdn(ig,ilay)=fdn(ig,ilay+1)+edn(ig,ilay)-ddn(ig,ilay)
    41             thetad(ig,ilay)=( fdn(ig,ilay+1)*thetad(ig,ilay+1) + edn(ig,ilay)*theta(ig,ilay) ) / (fdn(ig,ilay)+ddn(ig,ilay))
     51            if ( 1 == 0 ) then
     52               thetad(ig,ilay)=( fdn(ig,ilay+1)*thetad(ig,ilay+1) + edn(ig,ilay)*theta(ig,ilay) ) / (fdn(ig,ilay)+ddn(ig,ilay))
     53            else
     54               www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
     55               thetad(ig,ilay)=www*thetad(ig,ilay+1) + (1.-www)*theta(ig,ilay)
     56            endif
    4257         endif
    4358      enddo
    4459   enddo
    4560
    46 
     61   !Boucle pour l'updraft
     62   do ilay=1,nlay,1
     63      do ig=1,ngrid
     64         if (ilay.le.lmax(ig) .and. lmax(ig)>1) then
     65            fup(ig,ilay+1)=fup(ig,ilay)+eup(ig,ilay)-dup(ig,ilay)
     66            if (ilay == 1 ) then
     67               thetau(ig,ilay)=theta(ig,ilay)
     68            else
     69               !thetau(ig,ilay)=( fup(ig,ilay)*thetau(ig,ilay-1) + eup(ig,ilay)*theta(ig,ilay) ) / (fup(ig,ilay+1)+dup(ig,ilay))
     70               !eup(ig,ilay)=fup(ig,ilay+1)-fup(ig,ilay)+dup(ig,ilay)
     71               !thetau(ig,ilay)=( fup(ig,ilay)*thetau(ig,ilay-1) + (fup(ig,ilay+1)-fup(ig,ilay)+dup(ig,ilay))*theta(ig,ilay) ) / (fup(ig,ilay+1)+dup(ig,ilay))
     72               www=fup(ig,ilay)/(fup(ig,ilay+1)+dup(ig,ilay))
     73               !1-www=(fup(ig,ilay+1)+dup(ig,ilay)-fup(ig,ilay))/(fup(ig,ilay+1)+dup(ig,ilay))
     74               thetau(ig,ilay)=www*thetau(ig,ilay-1)+(1.-www)*theta(ig,ilay)
     75            endif
     76         endif
     77      enddo
     78   enddo
     79   !Boucle pour calculer le flux up
     80   do ilay=2,nlay,1
     81     do ig=1,ngrid
     82       fthu(ig,ilay)=fup(ig,ilay)*thetau(ig,ilay-1)
     83       fthd(ig,ilay)=-fdn(ig,ilay)*thetad(ig,ilay)
     84       !!!!ATTENTION HYPOTHESE de FLUX COMPENSATOIRE DESCENDANT ET DONC comme schema amont on va chercher theta au dessus!!!!!
     85       fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*theta(ig,ilay)
     86       fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay)
     87     enddo
     88   enddo
     89   !Boucle pour calculer theta
     90   do ilay=1,nlay,1
     91     do ig=1,ngrid
     92       dtheta(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
     93!       theta(ig,ilay)=theta(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
     94     enddo
     95   enddo
     96   if (1==0) then
     97    do ilay=1,nlay,1
     98     do ig=1,ngrid
     99       theta(ig,ilay)=theta(ig,ilay) + (fup(ig,ilay)*thetau(ig,ilay-1)-fup(ig,ilay+1)*thetau(ig,ilay) + &
     100       & (fup(ig,ilay+1)+fdn(ig,ilay+1))*theta(ig,ilay+1) - (fup(ig,ilay)+fdn(ig,ilay))*theta(ig,ilay) + &
     101       & fdn(ig,ilay+1)*thetad(ig,ilay+1)-fdn(ig,ilay)*thetad(ig,ilay))*(ptimestep/masse(ig,ilay))
     102     enddo
     103    enddo
     104   endif
     105! Il reste a coder :
     106! d(rho theta)/dt = - d/dz(rho w'theta')
     107! d(theta)/dt = -1/rho * d/dz(rho w'theta')
     108! hydrostatique : dp - rho g dz
     109! dz = - dp /  (rho g)
     110! 1 / dz = - (rho g ) / dp
     111
     112! d(theta)/dt = -1/rho * d/dp(rho w'theta') * (- rho g )
     113! d(theta)/dt =  d/dp(rho w'theta') * g
     114
     115
     116! ----> calculer d/dz(w'theta')
     117! w'theta' = alpha_up*(w_up-w_bar)*(theta_up-theta_bar)
     118!          + alpha_down*(w_down-w_bar)*(theta_down-theta_bar)
     119!          + (1-alpha_up-alpha_down)*(w_env-w_bar)*(theta_env-theta_bar)
     120! avec w_bar=0=alpha_up*w_up+alpha_down*w_dn+(1-alpha_up-alpha_down)w_env
     121! -> w_env= - (alpha_up*w_up+alpha_down*w_dn)/(1-alpha_up-alpha_down)
     122!  et on a : theta_bar=alpha_up*theta_up + alpha_down*theta_dn + (1-alpha_up-alpha_down)theta_env
     123!
     124! w'theta' = alpha_up*w_up*(theta_up-theta_bar)
     125!          + alpha_down*w_down*(theta_down-theta_bar)
     126!          + (1-alpha_up-alpha_down)*w_env*(theta_env-theta_bar)
     127
     128! rho*w'theta' = rho*alpha_up*w_up*(theta_up-theta_bar)
     129!              + rho*alpha_down*w_down*(theta_down-theta_bar)
     130!              + rho*(1-alpha_up-alpha_down)*w_env*(theta_env-theta_bar)
     131!
     132
     133! rho*w'theta' = fup*(theta_up-theta_bar)
     134!              + fdn*(theta_down-theta_bar)
     135!              - rho*(alpha_up*w_up+alpha_down*w_dn)*(theta_env-theta_bar)
     136
     137!              - rho*(alpha_up*w_up+alpha_down*w_dn)*((theta_bar-alpha_up*theta_up-alpha_down*theta_down) /(1-alpha_up-alpha_down)- theta_bar)
     138!              - rho*(alpha_up*w_up+alpha_down*w_dn)*(theta_bar*(1/(1-alpha_up-alpha_down)-1) - (alpha_up*theta_up+alpha_down*theta_down) /(1-alpha_up-alpha_down))
     139
     140! rho*w'theta' = fup*(theta_up-theta_bar)
     141!              + fdn*(theta_down-theta_bar)
     142!              - (fup+fdn)*(theta_env-theta_bar)
     143
     144
     145! rho*w'theta' = fup*(theta_up-theta_bar)
     146!              + fdn*(theta_down-theta_bar)
     147!              - (fup+fdn)*( (theta_bar-alpha_up*theta_up-alpha_down*theta_down) /(1-alpha_up-alpha_down) - theta_bar)
     148
     149!!! hypoth??se : alpha_up+alpha_down << 1 -> 1/(1-alpha_up-alpha_down) ~ 1
     150
     151! rho*w'theta' = fup*(theta_up-theta_bar)
     152
     153!              + fdn*(theta_down-theta_bar)
     154!              + (fup+fdn)*(alpha_up*theta_up+alpha_down*theta_down)
     155
     156! d(theta)/dt= -1/rho d/dz(rho*w'theta')
     157! choix de schema temporel (euler explicite)
     158! (theta(t,k)-theta(t-1,k))/dt = f(t-1,k)
     159! (theta(t,k)-theta(t-1,k))/dt = f(t,k) ---> euler implicite
     160! -----> on choisit explicite en temps
     161
     162!dans le cas d'une ascendance sans downdraft :
     163! d/dz(rho*w'theta') = fup(k)*qa(k-1) -fup(k+1)*qa(k)
     164
     165!(rho*w'theta'(k+1)-rho*w'theta'(k))/dz = fup(k)thetau(k-1)-fup(k+1)thetaup(k) + (fup+fdn)(k+1)*theta_env(k+1) - (fup+fdn)(k)*theta_env(k) + fdn(k+1) thetadn(k+1)
     166!                   - fdn(k) thetadn(k)
     167
     168! en continue on a  : d(theta)/dt = -1/rho * d/dz(rho w'theta')
     169! en discretis?? on a :
     170! (theta(t,k)-theta(t-1,k))/dt = -1/rho * rho*w'theta'(k+1)-rho*w'theta'(k))/dz
     171! theta(t,k) = theta(t-1,k) - dt/rho * (rho*w'theta'(k+1)-rho*w'theta'(k))/dz
     172! theta(t,k) = theta(t-1,k) - dt/rho * [fup(k)thetau(k-1)-fup(k+1)thetaup(k) + (fup+fdn)(k+1)*theta_env(k+1) -
     173! (fup+fdn)(k)*theta_env(k) + fdn(k+1) thetadn(k+1) - fdn(k) thetadn(k)]
     174
     175
     176
     177! d(theta)/dt =  d/dp(rho w'theta') * g
     178! -1/rho * d/dz(rho w'theta') = g*d/dp(rho w'theta') = g * [fup(k)thetau(k-1)-fup(k+1)thetaup(k) + (fup+fdn)(k+1)*theta_env(k+1) -
     179! (fup+fdn)(k)*theta_env(k) + fdn(k+1) thetadn(k+1) - fdn(k) thetadn(k)]
     180
     181
     182
     183! choix de sch??ma spatial
     184! d/dz(rho*w'theta') = (rho*w'theta'(k+1)-rho*w'theta'(k))*dz
     185! d/dz(rho*w'theta') = (rho*w'theta'(k)-rho*w'theta'(k-1))*dz
     186
     187
     188
     189! d/dz(rho*w'theta') = rho w'theta'(k+1)-rho w'theta'(k)
    47190
    48191!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    … …  
    61204
    62205!--------------------------------------------------------------
    63 !thermcell_down: calcul des propriétés du panache descendant.
     206!thermcell_down: calcul des propri??t??s du panache descendant.
    64207!--------------------------------------------------------------
    65208
    … …  
    98241! FH MODIFS APRES REUNIONS POUR COMMISSIONS
    99242! quelques erreurs de declaration
    100 ! probleme si lmax=1 ce qui a l'air d'être le cas en début de simu. Devrait être 0 ?
     243! probleme si lmax=1 ce qui a l'air d'??tre le cas en d??but de simu. Devrait ??tre 0 ?
    101244! Remarques :
    102 ! on pourrait écrire la formule de thetad
     245! on pourrait ??crire la formule de thetad
    103246!    www=fdn(ig,ilay+1)/ (fdn(ig,ilay)+ddn(ig,ilay))
    104247!    thetad(ig,ilay)= www * thetad(ig,ilay+1) + (1.-www) * theta(ig,ilay)
    105 ! Elle a l'avantage de bien montré la conservation, l'idée fondamentale dans le
     248! Elle a l'avantage de bien montr?? la conservation, l'id??e fondamentale dans le
    106249!   transport qu'on ne fait que sommer des "sources" au travers d'un "propagateur"
    107250!   (Green)
    108 ! Elle montre aussi beaucoup plus clairement pourquoi on n'a pas à se souccier (trop)
    109 !   de la possible nulité du dénominateur
     251! Elle montre aussi beaucoup plus clairement pourquoi on n'a pas ?? se souccier (trop)
     252!   de la possible nulit?? du d??nominateur
    110253
    111254
    … …  
    124267   ! Ecrire la conservervation de theta_l dans le panache descendant
    125268   ! Eventuellement faire la transformation theta_l -> theta_v
    126    ! Si l'air est sec (et qu'on oublie le côté theta_v) on peut
     269   ! Si l'air est sec (et qu'on oublie le c??t?? theta_v) on peut
    127270   ! se contenter de conserver theta.
    128271   !
    129    ! Connaissant thetadn, on peut calculer la flotabilité.
    130    ! Connaissant la flotabilité, on peut calculer w de proche en proche
    131    ! On peut calculer le detrainement de facon à garder alpha*rho = cste
    132    ! On en déduit l'entrainement latéral
    133    ! C'est le modèle des mini-projets.
     272   ! Connaissant thetadn, on peut calculer la flotabilit??.
     273   ! Connaissant la flotabilit??, on peut calculer w de proche en proche
     274   ! On peut calculer le detrainement de facon ?? garder alpha*rho = cste
     275   ! On en d??duit l'entrainement lat??ral
     276   ! C'est le mod??le des mini-projets.
    134277
    135278!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

  • LMDZ6/trunk/libf/phylmd/thermcell_main.F90

    r4365 r4377  
    1 
    21! $Id$
    32!
    … …  
    456455     &       detr,zqla,lev_out,lunout1,igout)
    457456
    458 #undef DevThermcellDown
    459 #ifdef DevThermcellDown
    460       print*,'WARNING !!! routine thermcell_down en cours de developpement'
    461       CALL thermcell_down(ngrid,nlay,po,pt,pu,pv,pplay,pplev,  &
    462      &           lmax,fm,entr,detr,zthl)
    463 
    464 #endif
    465 
    466457!IM 060508    &       detr,zqla,zmax,lev_out,lunout,igout)
    467458
    … …  
    490481!   calcul du transport vertical
    491482!------------------------------------------------------------------
     483
     484#undef DevThermcellDown
     485#ifdef DevThermcellDown
     486      print*,'AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'
     487      print*,'WARNING !!! routine thermcell_down en cours de developpement'
     488      CALL thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,entr0,detr0,0.5*detr0,0.5*entr0,masse,zthl)
     489!--------------------------------------------------------------
     490#endif
    492491
    493492      call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,  &

  • LMDZ6/trunk/makelmdz

    r4283 r4377  
    708708   # On enleve tout apres ONLy et on met un "uniq" pour que ca ne recrée pas
    709709   # le makefile si on se contente d'ajouter des lignes dans le ONLY
     710   exclude="replay automatic include"
    710711   for str in subroutine "use " "include " ; do
    711       grep -i "$str" libf/$dir/*.[Fh] | cut -d\( -f1 | sed -e 's/[Oo][Nn][Ll][Yy].*.$//' | uniq >> tmp77
    712       grep -i "$str" libf/$dir/*.F90  | cut -d\( -f1 | sed -e 's/[Oo][Nn][Ll][Yy].*.$//' | uniq >> tmp90
     712      grep -i "$str" libf/$dir/*.[Fh] | sed -e "/$exclude/d" | cut -d\( -f1 | sed -e 's/[Oo][Nn][Ll][Yy].*.$//' | uniq >> tmp77
     713      grep -i "$str" libf/$dir/*.F90  | sed -e "/$exclude/d" | cut -d\( -f1 | sed -e 's/[Oo][Nn][Ll][Yy].*.$//' | uniq >> tmp90
    713714   done
    714715done
Note: See TracChangeset for help on using the changeset viewer.