Ignore:
Timestamp:
Apr 29, 2019, 10:07:47 AM (6 years ago)
Author:
aboissinot
Message:

Fix in convadj.F
New version of the thermal plume model (new entrainment and detrainment formulae, alimentation is removed)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90

    r2102 r2127  
    22!
    33!
    4 SUBROUTINE thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,detr,masse,         &
    5                         q,dq,qa,lmin,lev_out)
     4SUBROUTINE thermcell_dq(ngrid,nlay,ptimestep,fm,entr,detr,masse,              &
     5                        q,dq,qa,lmin)
    66     
    77     
    8 !==============================================================================
    9 !  Calcul du transport verticale dans la couche limite en presence
    10 !  de "thermiques" explicitement representes
    11 !  calcul du dq/dt une fois qu'on connait les ascendances
    12 !
     8!===============================================================================
     9!  Purpose: Calcul du transport verticale dans la couche limite en presence de
     10!          "thermiques" explicitement representes
     11!           Calcul du dq/dt une fois qu'on connait les ascendances
     12! 
    1313!  Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
    14 !  Introduction of an implicit computation of vertical advection in
    15 !  the environment of thermal plumes in thermcell_dq
    16 !  impl =     0 : explicit, 1 : implicit, -1 : old version
    17 !
    18 !  Modif 2018/11/05 (AB alexandre.boissinot@lmd.jussieu.fr)
     14!  Introduction of an implicit computation of vertical advection in the environ-
     15!     ment of thermal plumes in thermcell_dq
    1916! 
    20 !
    21 !==============================================================================
     17!  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
     18!     dqimpl = 1 : implicit scheme
     19!     dqimpl = 0 : explicit scheme
     20! 
     21!===============================================================================
    2222     
    2323      USE print_control_mod, ONLY: prt_level
     24      USe thermcell_mod, ONLY: dqimpl
    2425     
    2526      IMPLICIT NONE
    2627     
    2728     
    28 !==============================================================================
     29!===============================================================================
    2930! Declaration
    30 !==============================================================================
     31!===============================================================================
    3132     
    32 !      inputs:
    33 !      -------
     33!     Inputs:
     34!     -------
    3435     
    3536      INTEGER ngrid, nlay
    36       INTEGER impl
    3737      INTEGER lmin(ngrid)
    38       INTEGER lev_out                           ! niveau pour les print
    3938     
    4039      REAL ptimestep
    … …  
    4443      REAL detr(ngrid,nlay)
    4544      REAL q(ngrid,nlay)
     45     
     46!     Outputs:
     47!     --------
     48     
     49      REAL dq(ngrid,nlay)
    4650      REAL qa(ngrid,nlay)
    4751     
    48 !      outputs:
    49 !      --------
    50      
    51       REAL dq(ngrid,nlay)
    52      
    53 !      local:
    54 !      ------
     52!     Local:
     53!     ------
    5554     
    5655      INTEGER ig, l
    … …  
    6362      REAL zzm
    6463     
    65       CHARACTER (LEN=20) :: modname
    66       CHARACTER (LEN=80) :: abort_message
     64!===============================================================================
     65! Initialization
     66!===============================================================================
    6767     
    68 !==============================================================================
    69 ! Initialization
    70 !==============================================================================
     68      qold(:,:) = q(:,:)
    7169     
    72       qold = q
     70!===============================================================================
     71! Tracer variation computation
     72!===============================================================================
    7373     
    74       modname = 'thermcell_dq'
    75      
    76       IF (impl.lt.0) THEN
    77          print *, 'OLD EXPLICIT SCHEME IS USED'
    78          CALL thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,               &
    79          &                   masse,q,dq,qa,lev_out)
    80          RETURN
    81       ENDIF
    82      
    83 !------------------------------------------------------------------------------
     74!-------------------------------------------------------------------------------
    8475! CFL criterion computation for advection in downdraft
    85 !------------------------------------------------------------------------------
     76!-------------------------------------------------------------------------------
    8677     
    8778      cfl = 0.
    … …  
    10394               print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)
    10495               print *, 'fm-', fm(ig,l-1)
    105                abort_message = 'entr dt > m, 1st'
    106                CALL abort_physic(modname,abort_message,1)
     96               CALL abort
    10797            ENDIF
    10898         ENDDO
    10999      ENDDO
    110100     
    111 !------------------------------------------------------------------------------
     101!-------------------------------------------------------------------------------
    112102! Computation of tracer concentrations in the ascending plume
    113 !------------------------------------------------------------------------------
     103!-------------------------------------------------------------------------------
    114104     
    115105      DO ig=1,ngrid
    … …  
    121111      DO ig=1,ngrid
    122112         DO l=lmin(ig)+1,nlay
    123             if ((fm(ig,l+1)+detr(ig,l))*ptimestep.gt.1.e-5*masse(ig,l)) then
     113            IF ((fm(ig,l+1)+detr(ig,l))*ptimestep.gt.1.e-6*masse(ig,l)) THEN
    124114               qa(ig,l) = (fm(ig,l) * qa(ig,l-1) + entr(ig,l) * q(ig,l))      &
    125115               &        / (fm(ig,l+1) + detr(ig,l))
    126             else
     116            ELSE
    127117               qa(ig,l) = q(ig,l)
    128             endif
     118            ENDIF
    129119           
    130120!            IF (qa(ig,l).lt.0.) THEN
    … …  
    140130      ENDDO
    141131     
    142 !------------------------------------------------------------------------------
    143 ! Plume vertical flux of q
    144 !------------------------------------------------------------------------------
     132!-------------------------------------------------------------------------------
     133! Plume vertical flux of tracer
     134!-------------------------------------------------------------------------------
    145135     
    146136      DO l=2,nlay-1
    … …  
    151141      fqa(:,nlay) = 0.
    152142     
    153 !------------------------------------------------------------------------------
     143!-------------------------------------------------------------------------------
    154144! Trace species evolution
    155 !------------------------------------------------------------------------------
     145!-------------------------------------------------------------------------------
    156146     
    157       IF (impl==0) THEN
    158          do l=1,nlay-1
     147      IF (dqimpl==0) THEN
     148         DO l=1,nlay-1
    159149            q(:,l) = q(:,l) + (fqa(:,l) - fqa(:,l+1) - fm(:,l) * q(:,l)       &
    160150            &      + fm(:,l+1) * q(:,l+1)) * ptimestep / masse(:,l)
    161          enddo
    162       ELSE
    163          do l=nlay-1,1,-1
     151         ENDDO
     152      ELSEIF (dqimpl==1) THEN
     153         DO l=nlay-1,1,-1
    164154            q(:,l) = ( q(:,l) + ptimestep / masse(:,l)                        &
    165155            &      * ( fqa(:,l) - fqa(:,l+1) + fm(:,l+1) * q(:,l+1) ) )       &
    166156            &      / ( 1. + fm(:,l) * ptimestep / masse(:,l) )
    167157         ENDDO
     158      ELSE
     159         print *, 'ERROR: No corresponding scheme for mixing computations!'
     160         print *, '       dqimpl must be equal to 1, 0 or -1 but'
     161         print *, 'dqimpl =', dqimpl
     162         call abort
    168163      ENDIF
    169164     
    170 !==============================================================================
     165!===============================================================================
    171166! Tendencies
    172 !==============================================================================
     167!===============================================================================
    173168     
    174169      DO l=1,nlay
    … …  
    182177RETURN
    183178END
    184 
    185 
    186 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    187 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    188 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    189 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    190 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    191 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    192 
    193 
    194 Subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,masse,            &
    195                           q,dq,qa,lev_out)
    196      
    197      
    198 !==============================================================================
    199 !
    200 !   Calcul du transport verticale dans la couche limite en presence
    201 !   de "thermiques" explicitement representes
    202 !   calcul du dq/dt une fois qu'on connait les ascendances
    203 !
    204 !==============================================================================
    205      
    206       USE print_control_mod, ONLY: prt_level
    207      
    208       implicit none
    209      
    210      
    211 !==============================================================================
    212 ! Declaration
    213 !==============================================================================
    214      
    215       integer ngrid,nlay,impl
    216 
    217       real ptimestep
    218       real masse(ngrid,nlay),fm(ngrid,nlay+1)
    219       real entr(ngrid,nlay)
    220       real q(ngrid,nlay)
    221       real dq(ngrid,nlay)
    222       integer lev_out                           ! niveau pour les print
    223 
    224       real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
    225 
    226       real zzm
    227 
    228       integer ig,l
    229       real cfl
    230 
    231       real qold(ngrid,nlay)
    232       real ztimestep
    233       integer niter,iter
    234       CHARACTER (LEN=20) :: modname='thermcell_dq'
    235       CHARACTER (LEN=80) :: abort_message
    236 
    237 
    238 
    239 ! Calcul du critere CFL pour l'advection dans la subsidence
    240       cfl = 0.
    241       do l=1,nlay
    242          do ig=1,ngrid
    243             zzm=masse(ig,l)/ptimestep
    244             cfl=max(cfl,fm(ig,l)/zzm)
    245             if (entr(ig,l).gt.zzm) then
    246                print*,'entr*dt>m,2',l,entr(ig,l)*ptimestep,masse(ig,l)
    247                abort_message = 'entr dt > m, 2nd'
    248                CALL abort_physic (modname,abort_message,1)
    249             endif
    250          enddo
    251       enddo
    252      
    253 !IM 090508     print*,'CFL CFL CFL CFL ',cfl
    254      
    255 #undef CFL
    256 #ifdef CFL
    257 ! On subdivise le calcul en niter pas de temps.
    258       niter=int(cfl)+1
    259 #else
    260       niter=1
    261 #endif
    262      
    263       ztimestep=ptimestep/niter
    264       qold=q
    265      
    266       do iter=1,niter
    267          if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
    268          
    269 !   calcul du detrainement
    270          do l=1,nlay
    271             do ig=1,ngrid
    272                detr(ig,l)=fm(ig,l)-fm(ig,l+1)+entr(ig,l)
    273 !               print*,'Q2 DQ ',detr(ig,l),fm(ig,l),entr(ig,l)
    274 !test
    275                if (detr(ig,l).lt.0.) then
    276                   entr(ig,l)=entr(ig,l)-detr(ig,l)
    277                   detr(ig,l)=0.
    278 !                  print*,'detr2<0!!!','ig=',ig,'l=',l,'f=',fm(ig,l),
    279 !     s            'f+1=',fm(ig,l+1),'e=',entr(ig,l),'d=',detr(ig,l)
    280                endif
    281                
    282 !               if (fm(ig,l+1).lt.0.) then
    283 !                  print*,'fm2<0!!!'
    284 !               endif
    285                
    286 !               if (entr(ig,l).lt.0.) then
    287 !                  print*,'entr2<0!!!'
    288 !               endif
    289             enddo
    290          enddo
    291          
    292 ! calcul de la valeur dans les ascendances
    293         do ig=1,ngrid
    294             qa(ig,1)=q(ig,1)
    295          enddo
    296          
    297          do l=2,nlay
    298             do ig=1,ngrid
    299                if ((fm(ig,l+1)+detr(ig,l))*ztimestep.gt.1.e-5*masse(ig,l)) then
    300                   qa(ig,l) = (fm(ig,l)*qa(ig,l-1)+entr(ig,l)*q(ig,l))         &
    301                   &        / (fm(ig,l+1)+detr(ig,l))
    302                else
    303                   qa(ig,l)=q(ig,l)
    304                endif
    305                
    306                if (qa(ig,l).lt.0.) then
    307 !                  print*,'qa<0!!!'
    308                endif
    309                
    310                if (q(ig,l).lt.0.) then
    311 !                  print*,'q<0!!!'
    312                endif
    313             enddo
    314          enddo
    315          
    316 ! Calcul du flux subsident
    317          
    318          do l=2,nlay
    319             do ig=1,ngrid
    320 #undef centre
    321 #ifdef centre
    322                 wqd(ig,l)=fm(ig,l)*0.5*(q(ig,l-1)+q(ig,l))
    323 #else
    324 
    325 #define plusqueun
    326 #ifdef plusqueun
    327 ! Schema avec advection sur plus qu'une maille.
    328                zzm=masse(ig,l)/ztimestep
    329                if (fm(ig,l)>zzm) then
    330                  wqd(ig,l)=zzm*q(ig,l)+(fm(ig,l)-zzm)*q(ig,l+1)
    331                else
    332                   wqd(ig,l)=fm(ig,l)*q(ig,l)
    333                endif
    334 #else
    335                wqd(ig,l)=fm(ig,l)*q(ig,l)
    336 #endif
    337 #endif
    338 
    339                if (wqd(ig,l).lt.0.) then
    340 !                  print*,'wqd<0!!!'
    341                endif
    342             enddo
    343          enddo
    344          
    345          do ig=1,ngrid
    346             wqd(ig,1)=0.
    347             wqd(ig,nlay+1)=0.
    348          enddo
    349          
    350 ! Calcul des tendances
    351          do l=1,nlay
    352             do ig=1,ngrid
    353                q(ig,l)=q(ig,l)+(detr(ig,l)*qa(ig,l)-entr(ig,l)*q(ig,l)           &
    354                &        -wqd(ig,l)+wqd(ig,l+1))                                  &
    355                &        *ztimestep/masse(ig,l)
    356 !               if (dq(ig,l).lt.0.) then
    357 !                  print*,'dq<0!!!'
    358 !               endif
    359            enddo
    360          enddo
    361       enddo
    362      
    363 ! Calcul des tendances
    364       do l=1,nlay
    365          do ig=1,ngrid
    366             dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep
    367             q(ig,l) = qold(ig,l)
    368          enddo
    369       enddo
    370      
    371      
    372 return
    373 end
Note: See TracChangeset for help on using the changeset viewer.