Changeset 4437 for LMDZ6/trunk/libf


Ignore:
Timestamp:
Feb 14, 2023, 5:52:31 PM (20 months ago)
Author:
fhourdin
Message:

Atelier downdrafts. Version Maelle. Suite

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r4396 r4437  
    1    SUBROUTINE thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,eup,dup,edn,ddn,masse,trac,dtrac)
     1SUBROUTINE thermcell_updown_dq(ngrid,nlay,ptimestep,lmax,eup,dup,edn,ddn,masse,trac,dtrac)
     2
     3USE thermcell_ini_mod, ONLY: iflag_thermals_down
     4
    25
    36!-----------------------------------------------------------------
    … …  
    3235! Local
    3336
    34    real, dimension(ngrid,nlay+1) :: fup,fdn,fthu,fthd,fthe,fthtot
    35    real, dimension(ngrid,nlay) :: tracu,tracd
    36    real :: www
     37   real, dimension(ngrid,nlay+1) :: fup,fdn,fc,fthu,fthd,fthe,fthtot
     38   real, dimension(ngrid,nlay) :: tracu,tracd,traci,tracold
     39   real :: www, mstar_inv
    3740   integer ig,ilay
    38 
     41   real, dimension(ngrid,nlay):: s1,s2,num !coefficients pour la resolution implicite
     42   integer :: iflag_impl=1 ! 0 pour explicite, 1 pour implicite "classique", 2 pour implicite avec entrainement et detrainement
     43   
    3944   fdn(:,:)=0.
    4045   fup(:,:)=0.
     46   fc(:,:)=0.
    4147   fthu(:,:)=0.
    4248   fthd(:,:)=0.
    … …  
    4551   tracd(:,:)=0.
    4652   tracu(:,:)=0.
     53   traci(:,:)=trac(:,:)
     54   tracold(:,:)=trac(:,:)
     55   s1(:,:)=0.
     56   s2(:,:)=0.
     57   num(:,:)=1.
     58
     59   if ( iflag_thermals_down < 10 ) then
     60        stop 'thermcell_down_dq = 0 or >= 10'
     61   else
     62        iflag_impl=iflag_thermals_down-10
     63   endif
     64     
    4765
    4866   ! lmax : indice tel que fu(kmax+1)=0
    … …  
    6482         endif
    6583      enddo
    66    enddo
     84   enddo !Fin boucle sur l'updraft
     85   fdn(:,1)=0.
    6786
    6887   !Boucle pour l'updraft
    … …  
    83102         endif
    84103      enddo
    85    enddo
    86    !Boucle pour calculer le flux du traceur flux updraft, flux downdraft, flux compensatoire
    87    do ilay=2,nlay,1
     104      enddo !fin boucle sur le downdraft
     105
     106   ! Calcul des flux des traceurs dans les updraft et les downdrfat
     107   ! et du flux de masse compensateur
     108   ! en ilay=1 et nlay+1, fthu=0 et fthd=0
     109   fthu(:,1)=0.
     110   fthu(:,nlay+1)=0.
     111   fthd(:,1)=0.
     112   fthd(:,nlay+1)=0.
     113   fc(:,1)=0.
     114   fc(:,nlay+1)=0.
     115   do ilay=2,nlay,1 !boucle sur les interfaces
    88116     do ig=1,ngrid
    89117       fthu(ig,ilay)=fup(ig,ilay)*tracu(ig,ilay-1)
    90118       fthd(ig,ilay)=-fdn(ig,ilay)*tracd(ig,ilay)
    91        !!!!ATTENTION HYPOTHESE de FLUX COMPENSATOIRE DESCENDANT ET DONC comme schema amont on va chercher trac au dessus!!!!!
    92        !!!! si ce n'est pas le cas on stoppe le code !!!!
    93        if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then
    94           write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq'
    95           stop         
    96        endif
    97        fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay)
    98        !! si on voulait le prendre en compte on
    99        !fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay-1)
    100        fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay)
     119       fc(ig,ilay)=fup(ig,ilay)-fdn(ig,ilay)
    101120     enddo
    102121   enddo
    103    !Boucle pour calculer trac
    104    do ilay=1,nlay,1
    105      do ig=1,ngrid
    106        dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))*(1./masse(ig,ilay))
    107 !       trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
    108      enddo
    109    enddo
    110 ! Il reste a coder :
    111 ! d(rho trac)/dt = - d/dz(rho w'trac')
    112 ! d(trac)/dt = -1/rho * d/dz(rho w'trac')
    113 ! hydrostatique : dp - rho g dz
    114 ! dz = - dp /  (rho g)
    115 ! 1 / dz = - (rho g ) / dp
    116 
    117 ! d(trac)/dt = -1/rho * d/dp(rho w'trac') * (- rho g )
    118 ! d(trac)/dt =  d/dp(rho w'trac') * g
    119 
    120 
    121 ! ----> calculer d/dz(w'trac')
    122 ! w'trac' = alpha_up*(w_up-w_bar)*(trac_up-trac_bar)
    123 !          + alpha_down*(w_down-w_bar)*(trac_down-trac_bar)
    124 !          + (1-alpha_up-alpha_down)*(w_env-w_bar)*(trac_env-trac_bar)
    125 ! avec w_bar=0=alpha_up*w_up+alpha_down*w_dn+(1-alpha_up-alpha_down)w_env
    126 ! -> w_env= - (alpha_up*w_up+alpha_down*w_dn)/(1-alpha_up-alpha_down)
    127 !  et on a : trac_bar=alpha_up*trac_up + alpha_down*trac_dn + (1-alpha_up-alpha_down)trac_env
    128 !
    129 ! w'trac' = alpha_up*w_up*(trac_up-trac_bar)
    130 !          + alpha_down*w_down*(trac_down-trac_bar)
    131 !          + (1-alpha_up-alpha_down)*w_env*(trac_env-trac_bar)
    132 
    133 ! rho*w'trac' = rho*alpha_up*w_up*(trac_up-trac_bar)
    134 !              + rho*alpha_down*w_down*(trac_down-trac_bar)
    135 !              + rho*(1-alpha_up-alpha_down)*w_env*(trac_env-trac_bar)
    136 !
    137 
    138 ! rho*w'trac' = fup*(trac_up-trac_bar)
    139 !              + fdn*(trac_down-trac_bar)
    140 !              - rho*(alpha_up*w_up+alpha_down*w_dn)*(trac_env-trac_bar)
    141 
    142 !              - rho*(alpha_up*w_up+alpha_down*w_dn)*((trac_bar-alpha_up*trac_up-alpha_down*trac_down) /(1-alpha_up-alpha_down)- trac_bar)
    143 !              - rho*(alpha_up*w_up+alpha_down*w_dn)*(trac_bar*(1/(1-alpha_up-alpha_down)-1) - (alpha_up*trac_up+alpha_down*trac_down) /(1-alpha_up-alpha_down))
    144 
    145 ! rho*w'trac' = fup*(trac_up-trac_bar)
    146 !              + fdn*(trac_down-trac_bar)
    147 !              - (fup+fdn)*(trac_env-trac_bar)
    148 
    149 
    150 ! rho*w'trac' = fup*(trac_up-trac_bar)
    151 !              + fdn*(trac_down-trac_bar)
    152 !              - (fup+fdn)*( (trac_bar-alpha_up*trac_up-alpha_down*trac_down) /(1-alpha_up-alpha_down) - trac_bar)
    153 
    154 !!! hypoth??se : alpha_up+alpha_down << 1 -> 1/(1-alpha_up-alpha_down) ~ 1
    155 
    156 ! rho*w'trac' = fup*(trac_up-trac_bar)
    157 
    158 !              + fdn*(trac_down-trac_bar)
    159 !              + (fup+fdn)*(alpha_up*trac_up+alpha_down*trac_down)
    160 
    161 ! d(trac)/dt= -1/rho d/dz(rho*w'trac')
    162 ! choix de schema temporel (euler explicite)
    163 ! (trac(t,k)-trac(t-1,k))/dt = f(t-1,k)
    164 ! (trac(t,k)-trac(t-1,k))/dt = f(t,k) ---> euler implicite
    165 ! -----> on choisit explicite en temps
    166 
    167 !dans le cas d'une ascendance sans downdraft :
    168 ! d/dz(rho*w'trac') = fup(k)*qa(k-1) -fup(k+1)*qa(k)
    169 
    170 !(rho*w'trac'(k+1)-rho*w'trac'(k))/dz = fup(k)tracu(k-1)-fup(k+1)tracup(k) + (fup+fdn)(k+1)*trac_env(k+1) - (fup+fdn)(k)*trac_env(k) + fdn(k+1) tracdn(k+1)
    171 !                   - fdn(k) tracdn(k)
    172 
    173 ! en continue on a  : d(trac)/dt = -1/rho * d/dz(rho w'trac')
    174 ! en discretis?? on a :
    175 ! (trac(t,k)-trac(t-1,k))/dt = -1/rho * rho*w'trac'(k+1)-rho*w'trac'(k))/dz
    176 ! trac(t,k) = trac(t-1,k) - dt/rho * (rho*w'trac'(k+1)-rho*w'trac'(k))/dz
    177 ! trac(t,k) = trac(t-1,k) - dt/rho * [fup(k)tracu(k-1)-fup(k+1)tracup(k) + (fup+fdn)(k+1)*trac_env(k+1) -
    178 ! (fup+fdn)(k)*trac_env(k) + fdn(k+1) tracdn(k+1) - fdn(k) tracdn(k)]
    179 
    180 
    181 
    182 ! d(trac)/dt =  d/dp(rho w'trac') * g
    183 ! -1/rho * d/dz(rho w'trac') = g*d/dp(rho w'trac') = g * [fup(k)tracu(k-1)-fup(k+1)tracup(k) + (fup+fdn)(k+1)*trac_env(k+1) -
    184 ! (fup+fdn)(k)*trac_env(k) + fdn(k+1) tracdn(k+1) - fdn(k) tracdn(k)]
    185 
    186 
    187 
    188 ! choix de sch??ma spatial
    189 ! d/dz(rho*w'trac') = (rho*w'trac'(k+1)-rho*w'trac'(k))*dz
    190 ! d/dz(rho*w'trac') = (rho*w'trac'(k)-rho*w'trac'(k-1))*dz
    191 
    192 
    193 
    194 ! d/dz(rho*w'trac') = rho w'trac'(k+1)-rho w'trac'(k)
    195 
    196 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    197 ! Initialisations :
    198 !------------------
    199 
    200 
    201 !
     122   
     123
     124   !Boucle pour calculer le flux du traceur flux updraft, flux downdraft, flux compensatoire
     125   !Methode explicite :
     126   if(iflag_impl==0) then
     127     do ilay=2,nlay,1
     128       do ig=1,ngrid
     129         !!!!ATTENTION HYPOTHESE de FLUX COMPENSATOIRE DESCENDANT ET DONC comme schema amont on va chercher trac au dessus!!!!!
     130         !!!! tentative de prise en compte d'un flux compensatoire montant  !!!!
     131         if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then
     132            write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq'
     133            stop         
     134            !fthe(ig,ilay)=(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay-1)
     135         else
     136            fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay)
     137         endif
     138         !! si on voulait le prendre en compte on
     139         !fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay-1)
     140         fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay)
     141       enddo
     142     enddo
     143     !Boucle pour calculer trac
     144     do ilay=1,nlay
     145       do ig=1,ngrid
     146         dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))/masse(ig,ilay)
     147!         trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
     148       enddo
     149     enddo !fin du calculer de la tendance du traceur avec la methode explicite
     150
     151   !!! Reecriture du schéma explicite avec les notations du schéma implicite
     152   else if(iflag_impl==-1) then
     153     write(*,*) 'nouveau schéma explicite !!!'
     154     !!! Calcul de s1
     155     do ilay=1,nlay
     156       do ig=1,ngrid
     157         s1(ig,ilay)=fthu(ig,ilay)-fthu(ig,ilay+1)+fthd(ig,ilay)-fthd(ig,ilay+1)
     158         s2(ig,ilay)=s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1)
     159       enddo
     160     enddo
     161
     162     do ilay=2,nlay,1
     163       do ig=1,ngrid
     164         if (fup(ig,ilay)-fdn(ig,ilay) .lt. 0.) then
     165            write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq'
     166            stop         
     167         else
     168            fthe(ig,ilay)=-(fup(ig,ilay)-fdn(ig,ilay))*trac(ig,ilay)
     169         endif
     170         fthtot(ig,ilay)=fthu(ig,ilay)+fthd(ig,ilay)+fthe(ig,ilay)
     171       enddo
     172     enddo
     173     !Boucle pour calculer trac
     174     do ilay=1,nlay
     175       do ig=1,ngrid
     176         ! dtrac(ig,ilay)=(fthtot(ig,ilay)-fthtot(ig,ilay+1))/masse(ig,ilay)
     177         dtrac(ig,ilay)=(s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1))/masse(ig,ilay)
     178!         trac(ig,ilay)=trac(ig,ilay) + (fthtot(ig,ilay)-fthtot(ig,ilay+1))*(ptimestep/masse(ig,ilay))
     179!         trac(ig,ilay)=trac(ig,ilay) + (s1(ig,ilay)+fthe(ig,ilay)-fthe(ig,ilay+1))*(ptimestep/masse(ig,ilay))
     180       enddo
     181     enddo !fin du calculer de la tendance du traceur avec la methode explicite
     182
     183   else if (iflag_impl==1) then
     184     write(*,*) 'SCHEMA IMPLICITE EN COURS DE DEVELOPPEMENT !'
     185     do ilay=1,nlay
     186       do ig=1,ngrid
     187         s1(ig,ilay)=fthu(ig,ilay)-fthu(ig,ilay+1)+fthd(ig,ilay)-fthd(ig,ilay+1)
     188       enddo
     189     enddo
     190     
     191     !Boucle pour calculer traci = trac((t+dt)
     192     do ilay=nlay-1,1,-1
     193       do ig=1,ngrid
     194         if((fup(ig,ilay)-fdn(ig,ilay)) .lt. 0) then
     195            write(*,*) 'flux compensatoire montant, cas non traite par thermcell_updown_dq dans le cas d une resolution implicite, ilay : ', ilay
     196            stop 
     197         else
     198           mstar_inv=ptimestep/masse(ig,ilay)
     199           traci(ig,ilay)=((traci(ig,ilay+1)*fc(ig,ilay+1)+s1(ig,ilay))*mstar_inv+tracold(ig,ilay))/(1.+fc(ig,ilay)*mstar_inv)
     200         endif
     201       enddo
     202     enddo
     203     do ilay=1,nlay
     204       do ig=1,ngrid
     205         dtrac(ig,ilay)=(traci(ig,ilay)-tracold(ig,ilay))/ptimestep
     206       enddo
     207     enddo
     208
     209   else
     210      write(*,*) 'valeur de iflag_impl non prevue'
     211      stop
     212
     213   endif
     214
    202215 RETURN
    203216   END
Note: See TracChangeset for help on using the changeset viewer.