Ignore:
Timestamp:
Feb 18, 2019, 11:40:47 AM (6 years ago)
Author:
aboissinot
Message:

Now detr is an argument of thermcell_dq subroutine (which is called in thermcell_main).
thermcell_dq is cleaned up.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2060 r2102  
    22!
    33!
    4       SUBROUTINE thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,              &
    5       &                       masse,q,dq,qa,lmin,lev_out)
    6      
    7      
    8       USE print_control_mod, ONLY: prt_level
    9      
    10       IMPLICIT NONE
     4SUBROUTINE thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,detr,masse,         &
     5                        q,dq,qa,lmin,lev_out)
     6     
    117     
    128!==============================================================================
    … …  
    2521!==============================================================================
    2622     
     23      USE print_control_mod, ONLY: prt_level
     24     
     25      IMPLICIT NONE
     26     
     27     
    2728!==============================================================================
    2829! Declaration
    … …  
    3233!      -------
    3334     
    34       INTEGER ngrid,nlay,impl
     35      INTEGER ngrid, nlay
     36      INTEGER impl
    3537      INTEGER lmin(ngrid)
    3638      INTEGER lev_out                           ! niveau pour les print
    3739     
    3840      REAL ptimestep
    39       REAL masse(ngrid,nlay),fm(ngrid,nlay+1)
     41      REAL masse(ngrid,nlay)
     42      REAL fm(ngrid,nlay+1)
    4043      REAL entr(ngrid,nlay)
     44      REAL detr(ngrid,nlay)
    4145      REAL q(ngrid,nlay)
    4246      REAL qa(ngrid,nlay)
    43       REAL detr(ngrid,nlay)
    4447     
    4548!      outputs:
    … …  
    5154!      ------
    5255     
    53       INTEGER ig,k
    54       INTEGER niter,iter
     56      INTEGER ig, l
     57      INTEGER niter, iter
    5558     
    5659      REAL cfl
    57       REAL qold(ngrid,nlay),fqa(ngrid,nlay+1)
     60      REAL qold(ngrid,nlay)
     61      REAL fqa(ngrid,nlay+1)
    5862      REAL wqd(ngrid,nlay+1)
    5963      REAL zzm
    … …  
    7882     
    7983!------------------------------------------------------------------------------
    80 ! Calcul du critere CFL pour l'advection dans la subsidence
     84! CFL criterion computation for advection in downdraft
    8185!------------------------------------------------------------------------------
    8286     
    8387      cfl = 0.
    8488     
    85       DO k=1,nlay
     89      DO l=1,nlay
    8690         DO ig=1,ngrid
    87             zzm = masse(ig,k) / ptimestep
    88             cfl = max(cfl, fm(ig,k) / zzm)
     91            zzm = masse(ig,l) / ptimestep
     92            cfl = max(cfl, fm(ig,l) / zzm)
    8993           
    90             IF (entr(ig,k).gt.zzm) THEN
     94            IF (entr(ig,l).gt.zzm) THEN
    9195               print *, 'ERROR: entrainment is greater than the layer mass!'
    92                print *, 'entr*dt,mass', entr(ig,k)*ptimestep, masse(ig,k)
    93                
    94                print *, 'ig,l', ig, k
    95                print *, 'fm-', fm(ig,k)
    96                print *, 'entr,detr', entr(ig,k), detr(ig,k)
    97                print *, 'fm+', fm(ig,k+1)
    98                
     96               print *, 'ig,l,entr', ig, l, entr(ig,l)
     97               print *, '-------------------------------'
     98               print *, 'entr*dt,mass', entr(ig,l)*ptimestep, masse(ig,l)
     99               print *, '-------------------------------'
     100               print *, 'fm+', fm(ig,l+1)
     101               print *, 'entr,detr', entr(ig,l), detr(ig,l)
     102               print *, 'fm ', fm(ig,l)
     103               print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)
     104               print *, 'fm-', fm(ig,l-1)
    99105               abort_message = 'entr dt > m, 1st'
    100106               CALL abort_physic(modname,abort_message,1)
    … …  
    104110     
    105111!------------------------------------------------------------------------------
    106 ! Detrainement computation
    107 !------------------------------------------------------------------------------
    108      
    109       DO k=1,nlay
    110          DO ig=1,ngrid
    111             detr(ig,k) = fm(ig,k) - fm(ig,k+1) + entr(ig,k)
     112! Computation of tracer concentrations in the ascending plume
     113!------------------------------------------------------------------------------
     114     
     115      DO ig=1,ngrid
     116         DO l=1,lmin(ig)
     117            qa(ig,l) = q(ig,l)
     118         ENDDO
     119      ENDDO
     120     
     121      DO ig=1,ngrid
     122         DO l=lmin(ig)+1,nlay
     123            if ((fm(ig,l+1)+detr(ig,l))*ptimestep.gt.1.e-5*masse(ig,l)) then
     124               qa(ig,l) = (fm(ig,l) * qa(ig,l-1) + entr(ig,l) * q(ig,l))      &
     125               &        / (fm(ig,l+1) + detr(ig,l))
     126            else
     127               qa(ig,l) = q(ig,l)
     128            endif
    112129           
    113             IF (detr(ig,k).lt.0.) THEN
    114                entr(ig,k) = entr(ig,k) - detr(ig,k)
    115                detr(ig,k) = 0.
    116 !               print *, 'WARNING: detr2 is negative!'
    117 !               print *, 'detr', detr(ig,k)
    118             ENDIF
    119            
    120 !            IF (qa(ig,k).lt.0.) THEN
    121 !               print *, 'WARNING: fm2 is negative!'
    122 !               print *, 'fm', fm(ig,k)
     130!            IF (qa(ig,l).lt.0.) THEN
     131!               print *, 'WARNING: qa is negative!'
     132!               print *, 'qa', qa(ig,l)
    123133!            ENDIF
    124134           
    125 !            IF (q(ig,k).lt.0.) THEN
    126 !               print *, 'WARNING: entr2 is negative!'
    127 !               print *, 'entr', entr(ig,k)
     135!            IF (q(ig,l).lt.0.) THEN
     136!               print *, 'WARNING: q is negative!'
     137!               print *, 'q', q(ig,l)
    128138!            ENDIF
    129139         ENDDO
    … …  
    131141     
    132142!------------------------------------------------------------------------------
    133 ! Computation of tracer concentrations in the ascending plume
    134 !------------------------------------------------------------------------------
    135      
    136       DO ig=1,ngrid
    137          DO k=1,lmin(ig)
    138             qa(ig,k) = q(ig,k)
    139          ENDDO
    140       ENDDO
    141      
    142       DO ig=1,ngrid
    143          DO k=lmin(ig)+1,nlay
    144             if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then
    145                qa(ig,k) = (fm(ig,k) * qa(ig,k-1) + entr(ig,k) * q(ig,k))      &
    146                &        / (fm(ig,k+1) + detr(ig,k))
    147             else
    148                qa(ig,k) = q(ig,k)
    149             endif
    150            
    151 !            IF (qa(ig,k).lt.0.) THEN
    152 !               print *, 'WARNING: qa is negative!'
    153 !               print *, 'qa', qa(ig,k)
    154 !            ENDIF
    155            
    156 !            IF (q(ig,k).lt.0.) THEN
    157 !               print *, 'WARNING: q is negative!'
    158 !               print *, 'q', q(ig,k)
    159 !            ENDIF
    160          ENDDO
    161       ENDDO
    162      
    163 !------------------------------------------------------------------------------
    164143! Plume vertical flux of q
    165144!------------------------------------------------------------------------------
    166145     
    167       DO k=2,nlay-1
    168          fqa(:,k) = fm(:,k) * qa(:,k-1)
     146      DO l=2,nlay-1
     147         fqa(:,l) = fm(:,l) * qa(:,l-1)
    169148      ENDDO
    170149     
    … …  
    177156     
    178157      IF (impl==0) THEN
    179          do k=1,nlay-1
    180             q(:,k) = q(:,k) + (fqa(:,k) - fqa(:,k+1) - fm(:,k) * q(:,k)       &
    181             &      + fm(:,k+1) * q(:,k+1)) * ptimestep / masse(:,k)
     158         do l=1,nlay-1
     159            q(:,l) = q(:,l) + (fqa(:,l) - fqa(:,l+1) - fm(:,l) * q(:,l)       &
     160            &      + fm(:,l+1) * q(:,l+1)) * ptimestep / masse(:,l)
    182161         enddo
    183162      ELSE
    184          do k=nlay-1,1,-1
    185             q(:,k) = ( q(:,k) + ptimestep / masse(:,k)                        &
    186             &      * ( fqa(:,k) - fqa(:,k+1) + fm(:,k+1) * q(:,k+1) ) )       &
    187             &      / ( 1. + fm(:,k) * ptimestep / masse(:,k) )
     163         do l=nlay-1,1,-1
     164            q(:,l) = ( q(:,l) + ptimestep / masse(:,l)                        &
     165            &      * ( fqa(:,l) - fqa(:,l+1) + fm(:,l+1) * q(:,l+1) ) )       &
     166            &      / ( 1. + fm(:,l) * ptimestep / masse(:,l) )
    188167         ENDDO
    189168      ENDIF
    … …  
    193172!==============================================================================
    194173     
    195       DO k=1,nlay
     174      DO l=1,nlay
    196175         DO ig=1,ngrid
    197             dq(ig,k) = (q(ig,k) - qold(ig,k)) / ptimestep
    198             q(ig,k) = qold(ig,k)
     176            dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep
     177            q(ig,l) = qold(ig,l)
    199178         ENDDO
    200179      ENDDO
    … …  
    205184
    206185
    207 
    208 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    209 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    210 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    211 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    212 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    213 
    214 
    215 
    216 ! Obsolete version kept for convergence with Cmip5 NPv3.1 simulations
    217 
    218 
    219       subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,  &
    220       &           masse,q,dq,qa,lev_out)
    221      
    222       USE print_control_mod, ONLY: prt_level
    223      
    224       implicit none
    225 
    226 !=======================================================================
     186!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     187!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     188!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     189!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     190!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     191!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     192
     193
     194Subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,masse,            &
     195                          q,dq,qa,lev_out)
     196     
     197     
     198!==============================================================================
    227199!
    228200!   Calcul du transport verticale dans la couche limite en presence
    … …  
    230202!   calcul du dq/dt une fois qu'on connait les ascendances
    231203!
    232 !=======================================================================
    233 
     204!==============================================================================
     205     
     206      USE print_control_mod, ONLY: prt_level
     207     
     208      implicit none
     209     
     210     
     211!==============================================================================
     212! Declaration
     213!==============================================================================
     214     
    234215      integer ngrid,nlay,impl
    235216
    … …  
    245226      real zzm
    246227
    247       integer ig,k
     228      integer ig,l
    248229      real cfl
    249230
    … …  
    258239! Calcul du critere CFL pour l'advection dans la subsidence
    259240      cfl = 0.
    260       do k=1,nlay
     241      do l=1,nlay
    261242         do ig=1,ngrid
    262             zzm=masse(ig,k)/ptimestep
    263             cfl=max(cfl,fm(ig,k)/zzm)
    264             if (entr(ig,k).gt.zzm) then
    265                print*,'entr*dt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k)
     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)
    266247               abort_message = 'entr dt > m, 2nd'
    267248               CALL abort_physic (modname,abort_message,1)
    … …  
    269250         enddo
    270251      enddo
    271 
     252     
    272253!IM 090508     print*,'CFL CFL CFL CFL ',cfl
    273 
     254     
    274255#undef CFL
    275256#ifdef CFL
    … …  
    279260      niter=1
    280261#endif
    281 
     262     
    282263      ztimestep=ptimestep/niter
    283264      qold=q
    284 
    285 
    286 do iter=1,niter
    287       if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
    288 
     265     
     266      do iter=1,niter
     267         if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
     268         
    289269!   calcul du detrainement
    290       do k=1,nlay
    291          do ig=1,ngrid
    292             detr(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
    293 !           print*,'Q2 DQ ',detr(ig,k),fm(ig,k),entr(ig,k)
     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)
    294274!test
    295             if (detr(ig,k).lt.0.) then
    296                entr(ig,k)=entr(ig,k)-detr(ig,k)
    297                detr(ig,k)=0.
    298 !               print*,'detr2<0!!!','ig=',ig,'k=',k,'f=',fm(ig,k),
    299 !     s         'f+1=',fm(ig,k+1),'e=',entr(ig,k),'d=',detr(ig,k)
    300             endif
    301             if (fm(ig,k+1).lt.0.) then
    302 !               print*,'fm2<0!!!'
    303             endif
    304             if (entr(ig,k).lt.0.) then
    305 !               print*,'entr2<0!!!'
    306             endif
    307          enddo
    308       enddo
    309 
    310 !   calcul de la valeur dans les ascendances
    311       do ig=1,ngrid
    312          qa(ig,1)=q(ig,1)
    313       enddo
    314 
    315       do k=2,nlay
    316          do ig=1,ngrid
    317             if ((fm(ig,k+1)+detr(ig,k))*ztimestep.gt.  &
    318      &         1.e-5*masse(ig,k)) then
    319          qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
    320      &         /(fm(ig,k+1)+detr(ig,k))
    321             else
    322                qa(ig,k)=q(ig,k)
    323             endif
    324             if (qa(ig,k).lt.0.) then
    325 !               print*,'qa<0!!!'
    326             endif
    327             if (q(ig,k).lt.0.) then
    328 !               print*,'q<0!!!'
    329             endif
    330          enddo
    331       enddo
    332 
     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         
    333316! Calcul du flux subsident
    334 
    335       do k=2,nlay
    336          do ig=1,ngrid
     317         
     318         do l=2,nlay
     319            do ig=1,ngrid
    337320#undef centre
    338321#ifdef centre
    339              wqd(ig,k)=fm(ig,k)*0.5*(q(ig,k-1)+q(ig,k))
     322                wqd(ig,l)=fm(ig,l)*0.5*(q(ig,l-1)+q(ig,l))
    340323#else
    341324
    … …  
    343326#ifdef plusqueun
    344327! Schema avec advection sur plus qu'une maille.
    345             zzm=masse(ig,k)/ztimestep
    346             if (fm(ig,k)>zzm) then
    347                wqd(ig,k)=zzm*q(ig,k)+(fm(ig,k)-zzm)*q(ig,k+1)
    348             else
    349                wqd(ig,k)=fm(ig,k)*q(ig,k)
    350             endif
     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
    351334#else
    352             wqd(ig,k)=fm(ig,k)*q(ig,k)
     335               wqd(ig,l)=fm(ig,l)*q(ig,l)
    353336#endif
    354337#endif
    355338
    356             if (wqd(ig,k).lt.0.) then
    357 !               print*,'wqd<0!!!'
    358             endif
     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
    359360         enddo
    360361      enddo
    361       do ig=1,ngrid
    362          wqd(ig,1)=0.
    363          wqd(ig,nlay+1)=0.
     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
    364369      enddo
    365      
    366 
    367 ! Calcul des tendances
    368       do k=1,nlay
    369          do ig=1,ngrid
    370             q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
    371      &               -wqd(ig,k)+wqd(ig,k+1))  &
    372      &               *ztimestep/masse(ig,k)
    373 !            if (dq(ig,k).lt.0.) then
    374 !               print*,'dq<0!!!'
    375 !            endif
    376          enddo
    377       enddo
    378 
    379 
    380 enddo
    381 
    382 
    383 ! Calcul des tendances
    384       do k=1,nlay
    385          do ig=1,ngrid
    386             dq(ig,k)=(q(ig,k)-qold(ig,k))/ptimestep
    387             q(ig,k)=qold(ig,k)
    388          enddo
    389       enddo
    390 
    391       return
    392       end
     370     
     371     
     372return
     373end
Note: See TracChangeset for help on using the changeset viewer.