Ignore:
Timestamp:
Feb 15, 2019, 2:43:57 PM (6 years ago)
Author:
aboissinot
Message:

Fix a bug in thermcell_alim.F90 where vertical and horizontal loops must be inverted.
In thermal plume model, arrays size declaration is standardised (no longer done thanks to dimphy module but by way of arguments).
Clean up some thermal plume model routines (remove uselesss variables...)

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2072 r2101  
    22!
    33!
    4       SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep                     &
    5       &                  ,pplay,pplev,pphi,debut                              &
    6       &                  ,pu,pv,pt,po                                         &
    7       &                  ,pduadj,pdvadj,pdtadj,pdoadj                         &
    8       &                  ,f0,fm0,entr0,detr0                                  &
    9       &                  ,zqta,zqla,ztv,ztva,ztla,zthl,zqsatth                &
    10       &                  ,zmax0,zw2,fraca                                     &
    11       &                  ,lmin,lmix,lalim,lmax                                &
    12       &                  ,zpspsk,ratqscth,ratqsdiff                           &
    13       &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th                    &
     4SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep                           &
     5                         ,pplay,pplev,pphi,firstcall                          &
     6                         ,pu,pv,pt,po                                         &
     7                         ,pduadj,pdvadj,pdtadj,pdoadj                         &
     8                         ,f0,fm0,entr0,detr0                                  &
     9                         ,zqta,zqla,ztv,ztva,ztla,zthl,zqsatth                &
     10                         ,zw2,fraca                                           &
     11                         ,lmin,lmix,lalim,lmax                                &
     12                         ,zpspsk,ratqscth,ratqsdiff                           &
     13                         ,Ale_bl,Alp_bl,lalim_conv,wght_th                    &
    1414!!! nrlmd le 10/04/2012
    15       &                  ,pbl_tke,pctsrf,omega,airephy                        &
    16       &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0   &
    17       &                  ,n2,s2,ale_bl_stat                                   &
    18       &                  ,therm_tke_max,env_tke_max                           &
    19       &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke          &
    20       &                  ,alp_bl_conv,alp_bl_stat)
     15                         ,pbl_tke,pctsrf,omega,airephy                        &
     16                         ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0   &
     17                         ,n2,s2,ale_bl_stat                                   &
     18                         ,therm_tke_max,env_tke_max                           &
     19                         ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke          &
     20                         ,alp_bl_conv,alp_bl_stat)
    2121!!! fin nrlmd le 10/04/2012
     22     
    2223     
    2324!==============================================================================
    … …  
    5152!==============================================================================
    5253     
    53       USE dimphy
    54 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    55 ! AB : why use dimphy to get klon, klev while subroutine arguments ngrid and
    56 !      nlay are the same?
    57 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    5854      USE thermcell_mod
    59       USE print_control_mod, ONLY: lunout,prt_level
    60       USE planete_mod, ONLY: preff
     55      USE print_control_mod, ONLY: lunout, prt_level
    6156     
    6257      IMPLICIT NONE
    … …  
    8479      REAL po(ngrid,nlay)                       !
    8580     
    86       LOGICAL debut
     81      LOGICAL firstcall
    8782     
    8883!   outputs:
    … …  
    9489      REAL pdoadj(ngrid,nlay)                   ! water convective variations
    9590     
    96       REAL fm0(klon,klev+1)                     ! mass flux      (after possible time relaxation)
    97       REAL entr0(klon,klev)                     ! entrainment    (after possible time relaxation)
    98       REAL detr0(klon,klev)                     ! detrainment    (after possible time relaxation)
    99       REAL f0(klon)                             ! mass flux norm (after possible time relaxation)
     91      REAL fm0(ngrid,nlay+1)                    ! mass flux      (after possible time relaxation)
     92      REAL entr0(ngrid,nlay)                    ! entrainment    (after possible time relaxation)
     93      REAL detr0(ngrid,nlay)                    ! detrainment    (after possible time relaxation)
     94      REAL f0(ngrid)                            ! mass flux norm (after possible time relaxation)
    10095     
    10196!   local:
    … …  
    118113     
    119114      INTEGER ig,k,l,ll,ierr
    120       INTEGER lmix_bis(klon)
    121       INTEGER lmax(klon)                        !
    122       INTEGER lmin(klon)                        !
    123       INTEGER lalim(klon)                       !
    124       INTEGER lmix(klon)                        !
    125      
    126       REAL linter(klon)
    127       REAL zmix(klon)
    128       REAL zmax(klon)
    129       REAL zw2(klon,klev+1)
    130       REAL zw_est(klon,klev+1)
    131       REAL zmax_sec(klon)
    132       REAL zmax0(klon)
    133      
    134       REAL zlay(klon,klev)                      ! layers altitude
    135       REAL zlev(klon,klev+1)                    ! levels altitude
    136       REAL rho(klon,klev)                       ! layers density
    137       REAL rhobarz(klon,klev)                   ! levels density
    138       REAL deltaz(klon,klev)                    ! layers heigth
    139       REAL masse(klon,klev)                     ! layers mass
    140       REAL zpspsk(klon,klev)                    ! Exner function
    141      
    142       REAL zu(klon,klev)                        ! environment zonal wind
    143       REAL zv(klon,klev)                        ! environment meridional wind
    144       REAL zo(klon,klev)                        ! environment water tracer
    145       REAL zva(klon,klev)                       ! plume zonal wind
    146       REAL zua(klon,klev)                       ! plume meridional wind
    147       REAL zoa(klon,klev)                       ! plume water tracer
    148      
    149       REAL zt(klon,klev)                        ! T    environment
    150       REAL zh(klon,klev)                        ! T,TP environment
    151       REAL zthl(klon,klev)                      ! TP   environment
    152       REAL ztv(klon,klev)                       ! TPV  environment ?
    153       REAL zl(klon,klev)                        ! ql   environment
    154      
    155       REAL zta(klon,klev)                       !
    156       REAL zha(klon,klev)                       ! TRPV plume
    157       REAL ztla(klon,klev)                      ! TP   plume
    158       REAL ztva(klon,klev)                      ! TRPV plume
    159       REAL ztva_est(klon,klev)                  ! TRPV plume (temporary)
    160       REAL zqla(klon,klev)                      ! qv   plume
    161       REAL zqta(klon,klev)                      ! qt   plume
    162      
    163       REAL wmax(klon)                           ! maximal vertical speed
    164       REAL wmax_tmp(klon)                       ! temporary maximal vertical speed
    165       REAL wmax_sec(klon)                       ! maximal vertical speed if dry case
    166      
    167       REAL fraca(klon,klev+1)                   ! updraft fraction
    168       REAL f_star(klon,klev+1)                  ! normalized mass flux
    169       REAL entr_star(klon,klev)                 ! normalized entrainment
    170       REAL detr_star(klon,klev)                 ! normalized detrainment
    171       REAL alim_star_tot(klon)                  ! integrated alimentation
    172       REAL alim_star(klon,klev)                 ! normalized alimentation
    173       REAL alim_star_clos(klon,klev)            ! closure alimentation
    174      
    175       REAL fm(klon,klev+1)                      ! mass flux
    176       REAL entr(klon,klev)                      ! entrainment
    177       REAL detr(klon,klev)                      ! detrainment
    178       REAL f(klon)                              ! mass flux norm
    179      
    180       REAL zdthladj(klon,klev)                  !
     115      INTEGER lmix_bis(ngrid)
     116      INTEGER lmax(ngrid)                       !
     117      INTEGER lmin(ngrid)                       !
     118      INTEGER lalim(ngrid)                      !
     119      INTEGER lmix(ngrid)                       !
     120     
     121      REAL linter(ngrid)
     122      REAL zmix(ngrid)
     123      REAL zmax(ngrid)
     124      REAL zw2(ngrid,nlay+1)
     125      REAL zw_est(ngrid,nlay+1)
     126      REAL zmax_sec(ngrid)
     127     
     128      REAL zlay(ngrid,nlay)                     ! layers altitude
     129      REAL zlev(ngrid,nlay+1)                   ! levels altitude
     130      REAL rho(ngrid,nlay)                      ! layers density
     131      REAL rhobarz(ngrid,nlay)                  ! levels density
     132      REAL deltaz(ngrid,nlay)                   ! layers heigth
     133      REAL masse(ngrid,nlay)                    ! layers mass
     134      REAL zpspsk(ngrid,nlay)                   ! Exner function
     135     
     136      REAL zu(ngrid,nlay)                       ! environment zonal wind
     137      REAL zv(ngrid,nlay)                       ! environment meridional wind
     138      REAL zo(ngrid,nlay)                       ! environment water tracer
     139      REAL zva(ngrid,nlay)                      ! plume zonal wind
     140      REAL zua(ngrid,nlay)                      ! plume meridional wind
     141      REAL zoa(ngrid,nlay)                      ! plume water tracer
     142     
     143      REAL zt(ngrid,nlay)                       ! T    environment
     144      REAL zh(ngrid,nlay)                       ! T,TP environment
     145      REAL zthl(ngrid,nlay)                     ! TP   environment
     146      REAL ztv(ngrid,nlay)                      ! TPV  environment ?
     147      REAL zl(ngrid,nlay)                       ! ql   environment
     148     
     149      REAL zta(ngrid,nlay)                      !
     150      REAL zha(ngrid,nlay)                      ! TRPV plume
     151      REAL ztla(ngrid,nlay)                     ! TP   plume
     152      REAL ztva(ngrid,nlay)                     ! TRPV plume
     153      REAL ztva_est(ngrid,nlay)                 ! TRPV plume (temporary)
     154      REAL zqla(ngrid,nlay)                     ! qv   plume
     155      REAL zqta(ngrid,nlay)                     ! qt   plume
     156     
     157      REAL wmax(ngrid)                          ! maximal vertical speed
     158      REAL wmax_tmp(ngrid)                      ! temporary maximal vertical speed
     159      REAL wmax_sec(ngrid)                      ! maximal vertical speed if dry case
     160     
     161      REAL fraca(ngrid,nlay+1)                  ! updraft fraction
     162      REAL f_star(ngrid,nlay+1)                 ! normalized mass flux
     163      REAL entr_star(ngrid,nlay)                ! normalized entrainment
     164      REAL detr_star(ngrid,nlay)                ! normalized detrainment
     165      REAL alim_star_tot(ngrid)                 ! integrated alimentation
     166      REAL alim_star(ngrid,nlay)                ! normalized alimentation
     167      REAL alim_star_clos(ngrid,nlay)           ! closure alimentation
     168     
     169      REAL fm(ngrid,nlay+1)                     ! mass flux
     170      REAL entr(ngrid,nlay)                     ! entrainment
     171      REAL detr(ngrid,nlay)                     ! detrainment
     172      REAL f(ngrid)                             ! mass flux norm
     173     
     174      REAL zdthladj(ngrid,nlay)                 !
    181175      REAL lambda                               ! time relaxation coefficent
    182176     
    183       REAL zsortie(klon,klev)
    184       REAL zsortie1d(klon)
     177      REAL zsortie(ngrid,nlay)
     178      REAL zsortie1d(ngrid)
    185179      REAL susqr2pi, Reuler
    186180      REAL zf
    187181      REAL zf2
    188       REAL thetath2(klon,klev)
    189       REAL wth2(klon,klev)
    190       REAL wth3(klon,klev)
    191       REAL q2(klon,klev)
     182      REAL thetath2(ngrid,nlay)
     183      REAL wth2(ngrid,nlay)
     184      REAL wth3(ngrid,nlay)
     185      REAL q2(ngrid,nlay)
    192186! FH : probleme de dimensionnement avec l'allocation dynamique
    193187!     common/comtherm/thetath2,wth2
    194       real wq(klon,klev)
    195       real wthl(klon,klev)
    196       real wthv(klon,klev)
    197       real ratqscth(klon,klev)
     188      real wq(ngrid,nlay)
     189      real wthl(ngrid,nlay)
     190      real wthv(ngrid,nlay)
     191      real ratqscth(ngrid,nlay)
    198192      real var
    199193      real vardiff
    200       real ratqsdiff(klon,klev)
     194      real ratqsdiff(ngrid,nlay)
    201195! niveau de condensation
    202       integer nivcon(klon)
    203       real zcon(klon)
     196      integer nivcon(ngrid)
     197      real zcon(ngrid)
    204198      real CHI
    205       real zcon2(klon)
    206       real pcon(klon)
    207       real zqsat(klon,klev)
    208       real zqsatth(klon,klev)
    209 ! FH/IM :   save f0
    210       real zlevinter(klon)
     199      real zcon2(ngrid)
     200      real pcon(ngrid)
     201      real zqsat(ngrid,nlay)
     202      real zqsatth(ngrid,nlay)
     203      real zlevinter(ngrid)
    211204      real seuil
    212205     
    … …  
    214207     
    215208!------Entrees
    216       real pbl_tke(klon,klev+1,nbsrf)
    217       real pctsrf(klon,nbsrf)
    218       real omega(klon,klev)
    219       real airephy(klon)
     209      real pbl_tke(ngrid,nlay+1,nbsrf)
     210      real pctsrf(ngrid,nbsrf)
     211      real omega(ngrid,nlay)
     212      real airephy(ngrid)
    220213!------Sorties
    221       real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
    222       real therm_tke_max0(klon),env_tke_max0(klon)
    223       real n2(klon),s2(klon)
    224       real ale_bl_stat(klon)
    225       real therm_tke_max(klon,klev),env_tke_max(klon,klev)
    226       real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
     214      real zlcl(ngrid),fraca0(ngrid),w0(ngrid),w_conv(ngrid)
     215      real therm_tke_max0(ngrid),env_tke_max0(ngrid)
     216      real n2(ngrid),s2(ngrid)
     217      real ale_bl_stat(ngrid)
     218      real therm_tke_max(ngrid,nlay),env_tke_max(ngrid,nlay)
     219      real alp_bl_det(ngrid),alp_bl_fluct_m(ngrid),alp_bl_fluct_tke(ngrid),alp_bl_conv(ngrid),alp_bl_stat(ngrid)
    227220!------Local
    228221      integer nsrf
    229       real rhobarz0(klon)                    ! Densite au LCL
    230       logical ok_lcl(klon)                   ! Existence du LCL des thermiques
    231       integer klcl(klon)                     ! Niveau du LCL
    232       real interp(klon)                      ! Coef d'interpolation pour le LCL
     222      real rhobarz0(ngrid)                      ! Densite au LCL
     223      logical ok_lcl(ngrid)                     ! Existence du LCL des thermiques
     224      integer klcl(ngrid)                       ! Niveau du LCL
     225      real interp(ngrid)                        ! Coef d'interpolation pour le LCL
    233226!------Triggering
    234       real Su                                ! Surface unite: celle d'un updraft elementaire
     227      real Su                                   ! Surface unite: celle d'un updraft elementaire
    235228      parameter(Su=4e4)
    236       real hcoef                             ! Coefficient directeur pour le calcul de s2
     229      real hcoef                                ! Coefficient directeur pour le calcul de s2
    237230      parameter(hcoef=1)
    238       real hmincoef                          ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2
     231      real hmincoef                             ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2
    239232      parameter(hmincoef=0.3)
    240       real eps1                              ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd)
     233      real eps1                                 ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd)
    241234      parameter(eps1=0.3)
    242       real hmin(ngrid)                       ! Ordonnee a l'origine pour le calcul de s2
    243       real zmax_moy(ngrid)                   ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
     235      real hmin(ngrid)                          ! Ordonnee a l'origine pour le calcul de s2
     236      real zmax_moy(ngrid)                      ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
    244237      real zmax_moy_coef
    245238      parameter(zmax_moy_coef=0.33)
    246       real depth(klon)                       ! Epaisseur moyenne du cumulus
    247       real w_max(klon)                       ! Vitesse max statistique
    248       real s_max(klon)
     239      real depth(ngrid)                         ! Epaisseur moyenne du cumulus
     240      real w_max(ngrid)                         ! Vitesse max statistique
     241      real s_max(ngrid)
    249242!------Closure
    250       real pbl_tke_max(klon,klev)            ! Profil de TKE moyenne
    251       real pbl_tke_max0(klon)                ! TKE moyenne au LCL
    252       real w_ls(klon,klev)                   ! Vitesse verticale grande echelle (m/s)
    253       real coef_m                            ! On considere un rendement pour alp_bl_fluct_m
     243      real pbl_tke_max(ngrid,nlay)              ! Profil de TKE moyenne
     244      real pbl_tke_max0(ngrid)                  ! TKE moyenne au LCL
     245      real w_ls(ngrid,nlay)                     ! Vitesse verticale grande echelle (m/s)
     246      real coef_m                               ! On considere un rendement pour alp_bl_fluct_m
    254247      parameter(coef_m=1.)
    255       real coef_tke                          ! On considere un rendement pour alp_bl_fluct_tke
     248      real coef_tke                             ! On considere un rendement pour alp_bl_fluct_tke
    256249      parameter(coef_tke=1.)
    257250     
    … …  
    259252     
    260253! Nouvelles variables pour la convection
    261       real Ale_bl(klon)
    262       real Alp_bl(klon)
    263       real alp_int(klon),dp_int(klon),zdp
    264       real ale_int(klon)
    265       integer n_int(klon)
    266       real fm_tot(klon)
    267       real wght_th(klon,klev)
    268       integer lalim_conv(klon)
     254      real Ale_bl(ngrid)
     255      real Alp_bl(ngrid)
     256      real alp_int(ngrid),dp_int(ngrid),zdp
     257      real ale_int(ngrid)
     258      integer n_int(ngrid)
     259      real fm_tot(ngrid)
     260      real wght_th(ngrid,nlay)
     261      integer lalim_conv(ngrid)
    269262     
    270263      CHARACTER*2 str2
    … …  
    282275      seuil = 0.25
    283276     
    284       IF (debut) THEN
     277      IF (firstcall) THEN
    285278         IF (iflag_thermals==15.or.iflag_thermals==16) THEN
    286279            dvimpl = 0
    … …  
    299292      entr(:,:) = 0.
    300293      detr(:,:) = 0.
    301      
    302       IF (ngrid.ne.klon) THEN
    303          print *, 'ERROR: ngrid and klon are different!'
    304          print *, 'ngrid =', ngrid
    305          print *, 'klon  =', klon
    306       ENDIF
    307      
    308       DO ig=1,klon
     294      f(:) = 0.
     295     
     296      DO ig=1,ngrid
    309297         f0(ig) = max(f0(ig), 1.e-2)
    310          zmax0(ig) = max(zmax0(ig),40.)
    311298      ENDDO
    312299     
    … …  
    352339     
    353340      zlev(:,1) = 0.
    354       zlev(:,nlay+1) = (2. * pphi(:,klev) - pphi(:,klev-1)) / RG
     341      zlev(:,nlay+1) = (2. * pphi(:,nlay) - pphi(:,nlay-1)) / RG
    355342     
    356343      DO l=1,nlay
    … …  
    467454!==============================================================================
    468455     
    469       CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,             &
    470       &    po,zl,rhobarz,zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,   &
    471       &    lalim,f0,detr_star,entr_star,f_star,ztva,                       &
    472       &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,   &
     456      CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,                &
     457      &    po,zl,rhobarz,zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,      &
     458      &    lalim,f0,detr_star,entr_star,f_star,ztva,                          &
     459      &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,      &
    473460      &    lmin,lev_out,lunout1,igout)
    474461     
    475       CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
    476       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
    477      
    478 !==============================================================================
    479 ! Calcul des caracteristiques du thermique:zmax,zmix,wmax
    480 !==============================================================================
    481      
    482       CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,   &
    483       &                              zlev,lmax,zmax,zmax0,zmix,wmax,f_star,   &
    484       &                              lev_out)
     462!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
     463!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
     464     
     465!==============================================================================
     466! Thermal plumes characteristics: zmax, zmix, wmax
     467!==============================================================================
     468     
     469      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,            &
     470      &                     zlev,lmax,zmax,zmix,wmax,f_star,lev_out)
    485471     
    486472!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    488474!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    489475     
    490       CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
    491       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
    492       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
    493       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
    494      
    495 !==============================================================================
    496 ! Closure and mass flux determining
     476!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
     477!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
     478!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
     479!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
     480     
     481!==============================================================================
     482! Closure and mass fluxes computation
    497483!==============================================================================
    498484     
    … …  
    500486      &                  lalim,lmin,zmax_sec,wmax_sec,lev_out)
    501487     
    502       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
    503       CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
     488!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
     489!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
    504490     
    505491!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    509495      alim_star_clos(:,:) = alim_star(:,:)
    510496      alim_star_clos(:,:) = entr_star(:,:) + alim_star(:,:)
    511          
     497     
    512498!------------------------------------------------------------------------------
    513499! Closure (dry if iflag_thermals_closure=1, moist if iflag_thermals_closure=2)
    … …  
    531517         f0 = (1.-lambda) * f + lambda * f0
    532518      ELSE
    533          f0 = f
     519         f0(:) = f(:)
    534520      ENDIF
    535521     
    … …  
    539525      IF (.not. (f0(1).ge.0.) ) THEN
    540526         abort_message = '.not. (f0(1).ge.0.)'
    541          write(*,*) 'f0 =', f0(1)
     527         print *, 'f0 =', f0(1)
    542528         CALL abort_physic(modname,abort_message,1)
     529      ELSE
     530         print *, 'f0 =', f0(1)
    543531      ENDIF
    544532     
    545 !==============================================================================
    546 ! Deduction des flux
    547 !==============================================================================
     533!------------------------------------------------------------------------------
     534! Mass fluxes
     535!------------------------------------------------------------------------------
    548536     
    549537      CALL thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
    550       &       lalim,lmin,lmax,alim_star,entr_star,detr_star,                  &
    551       &       f,rhobarz,zlev,zw2,fm,entr,                                     &
    552       &       detr,zqla,lev_out,lunout1,igout)
    553      
    554       CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
    555       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
     538      &                   lalim,lmin,lmax,alim_star,entr_star,detr_star,      &
     539      &                   f,rhobarz,zlev,zw2,fm,entr,detr,zqla,               &
     540      &                   lev_out,lunout1,igout)
     541     
     542!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
     543!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
    556544     
    557545!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    567555         detr0 = (1.-lambda) * detr + lambda * detr0
    568556      ELSE
    569          fm0   = fm
    570          entr0 = entr
    571          detr0 = detr
     557         fm0(:,:)   = fm(:,:)
     558         entr0(:,:) = entr(:,:)
     559         detr0(:,:) = detr(:,:)
    572560      ENDIF
    573561     
    574 !==============================================================================
    575 ! Calcul du transport vertical (de qt et tp)
    576 !==============================================================================
    577      
    578       CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
    579       &                 zthl,zdthladj,zta,lmin,lev_out)
    580      
    581       CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
    582       &                 po,pdoadj,zoa,lmin,lev_out)
    583      
    584       DO l=1,nlay
     562!------------------------------------------------------------------------------
     563! Updraft fraction
     564!------------------------------------------------------------------------------
     565     
     566      DO ig=1,ngrid
     567         fraca(ig,1) = 0.
     568         fraca(ig,nlay+1) = 0.
     569      ENDDO
     570     
     571      DO l=2,nlay
    585572         DO ig=1,ngrid
    586            pdtadj(ig,l) = zdthladj(ig,l) * zpspsk(ig,l) 
    587          ENDDO
    588       ENDDO
    589      
    590 !==============================================================================
    591 ! Calcul de la fraction de l'ascendance
    592 !==============================================================================
    593      
    594       DO ig=1,klon
    595          fraca(ig,1)=0.
    596          fraca(ig,nlay+1)=0.
    597       ENDDO
    598      
    599       DO l=2,nlay
    600          DO ig=1,klon
    601             IF (zw2(ig,l).gt.1.e-10) THEN
     573            IF (zw2(ig,l).gt.0.) THEN
    602574               fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l))
    603575            ELSE
    … …  
    608580     
    609581!==============================================================================
     582! Transport vertical
     583!==============================================================================
     584     
     585!------------------------------------------------------------------------------
     586! Calcul du transport vertical (de qt et tp)
     587!------------------------------------------------------------------------------
     588     
     589      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
     590      &                 zthl,zdthladj,zta,lmin,lev_out)
     591     
     592      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
     593      &                 po,pdoadj,zoa,lmin,lev_out)
     594     
     595      DO l=1,nlay
     596         DO ig=1,ngrid
     597           pdtadj(ig,l) = zdthladj(ig,l) * zpspsk(ig,l) 
     598         ENDDO
     599      ENDDO
     600     
     601!------------------------------------------------------------------------------
    610602! Calcul du transport vertical du moment horizontal
    611 !==============================================================================
     603!------------------------------------------------------------------------------
    612604     
    613605      IF (dvimpl.eq.0) THEN
    … …  
    620612      ELSE
    621613!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    622 ! Calcul purement conservatif pour le transport de V
     614! Calcul purement conservatif pour le transport de V=(u,v)
    623615!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    624616         CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,       &
    … …  
    756748         ratqsdiff(:,:) = 0.
    757749         
    758          DO l=1,klev
     750         DO l=1,nlay
    759751            DO ig=1,ngrid
    760752               IF (l<=lalim(ig)) THEN
    … …  
    766758         if (prt_level.ge.1) print*,'14f OK convect8'
    767759     
    768          DO l=1,klev
     760         DO l=1,nlay
    769761            DO ig=1,ngrid
    770762               IF (l<=lalim(ig)) THEN
    … …  
    799791
    800792
    801       SUBROUTINE test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,         &
    802                              ztva,zqla,f_star,zw2,comment)
     793SUBROUTINE test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po,              &
     794                       ztva,zqla,f_star,zw2,comment)
    803795     
    804796     
    … …  
    815807!      -------
    816808     
    817       INTEGER klon
    818       INTEGER klev
    819       INTEGER long(klon)
    820      
    821       REAL pplev(klon,klev+1)
    822       REAL pplay(klon,klev)
    823       REAL ztv(klon,klev)
    824       REAL po(klon,klev)
    825       REAL ztva(klon,klev)
    826       REAL zqla(klon,klev)
    827       REAL f_star(klon,klev)
    828       REAL zw2(klon,klev)
     809      INTEGER ngrid
     810      INTEGER nlay
     811      INTEGER long(ngrid)
     812     
     813      REAL pplev(ngrid,nlay+1)
     814      REAL pplay(ngrid,nlay)
     815      REAL ztv(ngrid,nlay)
     816      REAL po(ngrid,nlay)
     817      REAL ztva(ngrid,nlay)
     818      REAL zqla(ngrid,nlay)
     819      REAL f_star(ngrid,nlay)
     820      REAL zw2(ngrid,nlay)
    829821      REAL seuil
    830822     
    … …  
    847839     
    848840!  test sur la hauteur des thermiques ...
    849       do i=1,klon
     841      do i=1,ngrid
    850842!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
    851843        if (prt_level.ge.10) then
    852844            print *, 'WARNING ',comment,' au point ',i,' K= ',long(i)
    853845            print *, '  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
    854             do k=1,klev
     846            do k=1,nlay
    855847               write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
    856848            enddo
    … …  
    871863
    872864
    873       SUBROUTINE thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,      &
    874       &                                  rg,pplev,therm_tke_max)
    875      
    876      
    877       USE print_control_mod, ONLY: prt_level
    878      
    879       IMPLICIT NONE
     865SUBROUTINE thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,            &
     866                                   rg,pplev,therm_tke_max)
    880867     
    881868     
    … …  
    887874!   Transport de la TKE par le thermique moyen pour la fermeture en ALP
    888875!   On transporte pbl_tke pour donner therm_tke
     876!==============================================================================
     877     
     878      USE print_control_mod, ONLY: prt_level
     879     
     880      IMPLICIT NONE
     881     
     882     
     883!==============================================================================
     884! Declarations
    889885!==============================================================================
    890886     
    … …  
    900896      real entr(ngrid,nlay)
    901897      real q(ngrid,nlay)
    902       integer lev_out                           ! niveau pour les print
    903898     
    904899      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
    … …  
    909904      integer isrf
    910905     
    911       lev_out=0
    912      
    913       if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
    914      
    915 !   calcul du detrainement
     906!------------------------------------------------------------------------------
     907! Calcul du detrainement
     908!------------------------------------------------------------------------------
     909     
    916910      do k=1,nlay
    917          detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
    918          masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
     911         detr0(:,k) = fm0(:,k) - fm0(:,k+1) + entr0(:,k)
     912         masse0(:,k) = (pplev(:,k) - pplev(:,k+1)) / RG
    919913      enddo
    920 
    921 
     914     
     915!------------------------------------------------------------------------------
    922916! Decalage vertical des entrainements et detrainements.
     917!------------------------------------------------------------------------------
     918     
    923919      masse(:,1)=0.5*masse0(:,1)
    924920      entr(:,1)=0.5*entr0(:,1)
    925921      detr(:,1)=0.5*detr0(:,1)
    926922      fm(:,1)=0.
     923     
    927924      do k=1,nlay-1
    928925         masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
    … …  
    931928         fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
    932929      enddo
     930     
    933931      fm(:,nlay+1)=0.
    934 
     932     
    935933!!! nrlmd le 16/09/2010
    936934!   calcul de la valeur dans les ascendances
    937 !       do ig=1,ngrid
    938 !          qa(ig,1)=q(ig,1)
    939 !       enddo
    940 !!!
    941 
    942 !do isrf=1,nsrf
    943 
    944 !   q(:,:)=therm_tke(:,:,isrf)
    945    q(:,:)=therm_tke_max(:,:)
    946 !!! nrlmd le 16/09/2010
     935!      do ig=1,ngrid
     936!         qa(ig,1) = q(ig,1)
     937!      enddo
     938     
     939      q(:,:)=therm_tke_max(:,:)
     940     
    947941      do ig=1,ngrid
    948942         qa(ig,1)=q(ig,1)
    949943      enddo
    950 !!!
    951 
    952     if (1==1) then
    953       do k=2,nlay
     944     
     945      if (1==1) then
     946         do k=2,nlay
     947            do ig=1,ngrid
     948               if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then
     949                  qa(ig,k) = (fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))         &
     950                  &        / (fm(ig,k+1)+detr(ig,k))
     951               else
     952                  qa(ig,k)=q(ig,k)
     953               endif
     954               
     955!               if (qa(ig,k).lt.0.) print *, 'WARNING: qa is negative!'
     956!               if (q(ig,k).lt.0.) print *, 'WARNING: q is negative!'
     957            enddo
     958         enddo
     959         
     960!------------------------------------------------------------------------------
     961! Calcul du flux subsident
     962!------------------------------------------------------------------------------
     963         
     964         do k=2,nlay
     965            do ig=1,ngrid
     966               wqd(ig,k) = fm(ig,k) * q(ig,k)
     967               if (wqd(ig,k).lt.0.) print*,'WARNING: wqd is negative!'
     968            enddo
     969         enddo
     970         
    954971         do ig=1,ngrid
    955             if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
    956      &         1.e-5*masse(ig,k)) then
    957          qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
    958      &         /(fm(ig,k+1)+detr(ig,k))
    959             else
    960                qa(ig,k)=q(ig,k)
    961             endif
    962             if (qa(ig,k).lt.0.) then
    963 !               print*,'qa<0!!!'
    964             endif
    965             if (q(ig,k).lt.0.) then
    966 !               print*,'q<0!!!'
    967             endif
    968          enddo
    969       enddo
    970 
    971 ! Calcul du flux subsident
    972 
    973       do k=2,nlay
    974          do ig=1,ngrid
    975             wqd(ig,k)=fm(ig,k)*q(ig,k)
    976             if (wqd(ig,k).lt.0.) then
    977 !               print*,'wqd<0!!!'
    978             endif
    979          enddo
    980       enddo
    981       do ig=1,ngrid
    982          wqd(ig,1)=0.
    983          wqd(ig,nlay+1)=0.
    984       enddo
    985 
     972            wqd(ig,1) = 0.
     973            wqd(ig,nlay+1) = 0.
     974         enddo
     975         
     976!------------------------------------------------------------------------------
    986977! Calcul des tendances
    987       do k=1,nlay
    988          do ig=1,ngrid
    989             q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
    990      &               -wqd(ig,k)+wqd(ig,k+1))  &
    991      &               *ptimestep/masse(ig,k)
    992          enddo
    993       enddo
    994 
    995  endif
    996 
    997    therm_tke_max(:,:)=q(:,:)
     978!------------------------------------------------------------------------------
     979         
     980         do k=1,nlay
     981            do ig=1,ngrid
     982               q(ig,k) = q(ig,k) + ptimestep / masse(ig,k)                    &
     983               &       * (detr(ig,k) * qa(ig,k) - entr(ig,k) * q(ig,k)        &
     984               &       - wqd(ig,k) + wqd(ig,k+1))
     985            enddo
     986         enddo
     987      endif
     988     
     989      therm_tke_max(:,:) = q(:,:)
    998990     
    999991     
Note: See TracChangeset for help on using the changeset viewer.