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_main.F90

    r2101 r2102  
    22!
    33!
    4 SUBROUTINE 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                    &
     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
    2222     
    … …  
    213213!------Sorties
    214214      real zlcl(ngrid),fraca0(ngrid),w0(ngrid),w_conv(ngrid)
    215       real therm_tke_max0(ngrid),env_tke_max0(ngrid)
     215      real therm_tke_max0(ngrid)
     216      real env_tke_max0(ngrid)
    216217      real n2(ngrid),s2(ngrid)
    217218      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)
     219      real therm_tke_max(ngrid,nlay)
     220      real env_tke_max(ngrid,nlay)
     221      real alp_bl_det(ngrid)
     222      real alp_bl_fluct_m(ngrid)
     223      real alp_bl_fluct_tke(ngrid)
     224      real alp_bl_conv(ngrid)
     225      real alp_bl_stat(ngrid)
    220226!------Local
    221227      integer nsrf
    … …  
    225231      real interp(ngrid)                        ! Coef d'interpolation pour le LCL
    226232!------Triggering
    227       real Su                                   ! Surface unite: celle d'un updraft elementaire
    228       parameter(Su=4e4)
    229       real hcoef                                ! Coefficient directeur pour le calcul de s2
    230       parameter(hcoef=1)
    231       real hmincoef                             ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2
    232       parameter(hmincoef=0.3)
    233       real eps1                                 ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd)
    234       parameter(eps1=0.3)
     233      real,parameter :: Su=4e4                  ! Surface unite: celle d'un updraft elementaire
     234      real,parameter :: hcoef                   ! Coefficient directeur pour le calcul de s2
     235      real,parameter :: hmincoef=0.3            ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2 (hmincoef=0.3)
     236      real,parameter :: eps1=0.3                ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd)
    235237      real hmin(ngrid)                          ! Ordonnee a l'origine pour le calcul de s2
    236238      real zmax_moy(ngrid)                      ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
    237       real zmax_moy_coef
    238       parameter(zmax_moy_coef=0.33)
     239      real,parameter :: zmax_moy_coef=0.33
    239240      real depth(ngrid)                         ! Epaisseur moyenne du cumulus
    240241      real w_max(ngrid)                         ! Vitesse max statistique
    … …  
    244245      real pbl_tke_max0(ngrid)                  ! TKE moyenne au LCL
    245246      real w_ls(ngrid,nlay)                     ! Vitesse verticale grande echelle (m/s)
    246       real coef_m                               ! On considere un rendement pour alp_bl_fluct_m
    247       parameter(coef_m=1.)
    248       real coef_tke                             ! On considere un rendement pour alp_bl_fluct_tke
    249       parameter(coef_tke=1.)
     247      real,parameter :: coef_m=1.               ! On considere un rendement pour alp_bl_fluct_m
     248      real,parameter :: coef_tke=1.             ! On considere un rendement pour alp_bl_fluct_tke
    250249     
    251250!!! fin nrlmd le 10/04/2012
    252251     
    253252! Nouvelles variables pour la convection
     253      integer lalim_conv(ngrid)
     254      integer n_int(ngrid)
    254255      real Ale_bl(ngrid)
    255256      real Alp_bl(ngrid)
    256       real alp_int(ngrid),dp_int(ngrid),zdp
     257      real alp_int(ngrid)
     258      real dp_int(ngrid),zdp
    257259      real ale_int(ngrid)
    258       integer n_int(ngrid)
    259260      real fm_tot(ngrid)
    260261      real wght_th(ngrid,nlay)
    261       integer lalim_conv(ngrid)
    262262     
    263263      CHARACTER*2 str2
    … …  
    515515      IF (tau_thermals>1.) THEN
    516516         lambda = exp(-ptimestep/tau_thermals)
    517          f0 = (1.-lambda) * f + lambda * f0
     517         f0(:) = (1.-lambda) * f(:) + lambda * f0(:)
    518518      ELSE
    519519         f0(:) = f(:)
    … …  
    527527         print *, 'f0 =', f0(1)
    528528         CALL abort_physic(modname,abort_message,1)
    529       ELSE
    530          print *, 'f0 =', f0(1)
    531529      ENDIF
    532530     
    … …  
    587585!------------------------------------------------------------------------------
    588586     
    589       CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
     587      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse,    &
    590588      &                 zthl,zdthladj,zta,lmin,lev_out)
    591589     
    592       CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
     590      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse,    &
    593591      &                 po,pdoadj,zoa,lmin,lev_out)
    594592     
    … …  
    614612! Calcul purement conservatif pour le transport de V=(u,v)
    615613!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    616          CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,      &
     614         CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
    617615         &                 zu,pduadj,zua,lmin,lev_out)
    618616         
    619          CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,      &
     617         CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
    620618         &                 zv,pdvadj,zva,lmin,lev_out)
    621619      ENDIF
    … …  
    631629!------------------------------------------------------------------------------
    632630         
    633          if (prt_level.ge.1) print*,'14a OK convect8'
    634          
    635          do ig=1,ngrid
     631         DO ig=1,ngrid
    636632            nivcon(ig) = 0
    637633            zcon(ig) = 0.
    638          enddo
     634         ENDDO
    639635!nouveau calcul
    640636         do ig=1,ngrid
    641637            CHI = zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
    642638            pcon(ig) = pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
    643          enddo
     639         ENDDO
    644640         
    645641!IM       do k=1,nlay
    646642         do k=1,nlay-1
    647643            do ig=1,ngrid
    648                if ((pcon(ig).le.pplay(ig,k)).and.(pcon(ig).gt.pplay(ig,k+1))) then
    649                   zcon2(ig) = zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100.
    650                endif
    651             enddo
    652          enddo
     644               IF ((pcon(ig).le.pplay(ig,k)).and.(pcon(ig).gt.pplay(ig,k+1))) then
     645                  zcon2(ig) = zlay(ig,k) - (pcon(ig) - pplay(ig,k))           &
     646                            / (RG * rho(ig,k)) / 100.
     647               ENDIF
     648            ENDDO
     649         ENDDO
    653650         
    654651         ierr = 0
    655652         
    656653         do ig=1,ngrid
    657            if (pcon(ig).le.pplay(ig,nlay)) then
    658               zcon2(ig) = zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100.
     654           IF (pcon(ig).le.pplay(ig,nlay)) then
     655              zcon2(ig) = zlay(ig,nlay) - (pcon(ig) - pplay(ig,nlay))         &
     656                        / (RG * rho(ig,nlay)) / 100.
    659657              ierr = 1
    660             endif
    661          enddo
    662          
    663          if (ierr==1) then
     658            ENDIF
     659         ENDDO
     660         
     661         IF (ierr==1) then
    664662              write(*,*) 'ERROR: thermal plumes rise the model top'
    665663              CALL abort
    666          endif
    667          
    668          if (prt_level.ge.1) print*,'14b OK convect8'
     664         ENDIF
     665         
     666         IF (prt_level.ge.1) print*,'14b OK convect8'
    669667         
    670668         do k=nlay,1,-1
    671669            do ig=1,ngrid
    672                if (zqla(ig,k).gt.1e-10) then
     670               IF (zqla(ig,k).gt.1e-10) then
    673671                  nivcon(ig) = k
    674672                  zcon(ig) = zlev(ig,k)
    675                endif
    676             enddo
    677          enddo
    678          
    679          if (prt_level.ge.1) print*,'14c OK convect8'
     673               ENDIF
     674            ENDDO
     675         ENDDO
     676         
     677         IF (prt_level.ge.1) print*,'14c OK convect8'
    680678         
    681679!------------------------------------------------------------------------------
    … …  
    690688               ratqscth(ig,l) = 0.
    691689               ratqsdiff(ig,l) = 0.
    692             enddo
    693          enddo
    694          
    695          if (prt_level.ge.1) print*,'14d OK convect8'
     690            ENDDO
     691         ENDDO
     692         
     693         IF (prt_level.ge.1) print*,'14d OK convect8'
    696694         
    697695         do l=1,nlay
    … …  
    701699               thetath2(ig,l) = zf2*(ztla(ig,l)-zthl(ig,l))**2
    702700               
    703                if (zw2(ig,l).gt.1.e-10) then
     701               IF (zw2(ig,l).gt.1.e-10) then
    704702                  wth2(ig,l) = zf2*(zw2(ig,l))**2
    705703               else
    706704                  wth2(ig,l) = 0.
    707                endif
     705               ENDIF
    708706               
    709                wth3(ig,l) = zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))               &
     707               wth3(ig,l) = zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))            &
    710708               &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
    711709               q2(ig,l) = zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
    712710!test: on calcul q2/po=ratqsc
    713711               ratqscth(ig,l) = sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
    714             enddo
    715          enddo
     712            ENDDO
     713         ENDDO
    716714         
    717715!------------------------------------------------------------------------------
    … …  
    724722               wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
    725723               wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
    726             enddo
    727          enddo
    728          
    729          CALL thermcell_alp(ngrid,nlay,ptimestep,                                &
    730          &                  pplay,pplev,                                         &
    731          &                  fm0,entr0,lmax,                                      &
    732          &                  Ale_bl,Alp_bl,lalim_conv,wght_th,                    &
    733          &                  zw2,fraca,                                           &
    734          &                  pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,  &
    735          &                  pbl_tke,pctsrf,omega,airephy,                        &
    736          &                  zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,   &
    737          &                  n2,s2,ale_bl_stat,                                   &
    738          &                  therm_tke_max,env_tke_max,                           &
    739          &                  alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke,          &
     724            ENDDO
     725         ENDDO
     726         
     727         CALL thermcell_alp(ngrid,nlay,ptimestep,                             &
     728         &                  pplay,pplev,                                      &
     729         &                  fm0,entr0,lmax,                                   &
     730         &                  Ale_bl,Alp_bl,lalim_conv,wght_th,                 &
     731         &                  zw2,fraca,pcon,                                   &
     732         &                  rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,    &
     733         &                  pbl_tke,pctsrf,omega,airephy,                     &
     734         &                  zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,&
     735         &                  n2,s2,ale_bl_stat,                                &
     736         &                  therm_tke_max,env_tke_max,                        &
     737         &                  alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke,       &
    740738         &                  alp_bl_conv,alp_bl_stat)
    741739         
    … …  
    756754         ENDDO
    757755         
    758          if (prt_level.ge.1) print*,'14f OK convect8'
     756         IF (prt_level.ge.1) print*,'14f OK convect8'
    759757     
    760758         DO l=1,nlay
    … …  
    763761                  zf  = fraca(ig,l)
    764762                  zf2 = zf / (1.-zf)
    765                   vardiff = vardiff + alim_star(ig,l) * (zqta(ig,l) * 1000. - var)**2
     763                  vardiff = vardiff + alim_star(ig,l)                         &
     764                          * (zqta(ig,l) * 1000. - var)**2
    766765               ENDIF
    767766            ENDDO
    … …  
    840839!  test sur la hauteur des thermiques ...
    841840      do i=1,ngrid
    842 !IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
    843         if (prt_level.ge.10) then
     841!IMtemp IF (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
     842        IF (prt_level.ge.10) then
    844843            print *, 'WARNING ',comment,' au point ',i,' K= ',long(i)
    845844            print *, '  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
    846845            do k=1,nlay
    847846               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)
    848             enddo
    849         endif
    850       enddo
     847            ENDDO
     848        ENDIF
     849      ENDDO
    851850     
    852851     
    … …  
    911910         detr0(:,k) = fm0(:,k) - fm0(:,k+1) + entr0(:,k)
    912911         masse0(:,k) = (pplev(:,k) - pplev(:,k+1)) / RG
    913       enddo
     912      ENDDO
    914913     
    915914!------------------------------------------------------------------------------
    … …  
    927926         detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
    928927         fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
    929       enddo
     928      ENDDO
    930929     
    931930      fm(:,nlay+1)=0.
    … …  
    935934!      do ig=1,ngrid
    936935!         qa(ig,1) = q(ig,1)
    937 !      enddo
     936!      ENDDO
    938937     
    939938      q(:,:)=therm_tke_max(:,:)
    … …  
    941940      do ig=1,ngrid
    942941         qa(ig,1)=q(ig,1)
    943       enddo
    944      
    945       if (1==1) then
     942      ENDDO
     943     
     944      IF (1==1) then
    946945         do k=2,nlay
    947946            do ig=1,ngrid
    948                if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then
     947               IF ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then
    949948                  qa(ig,k) = (fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))         &
    950949                  &        / (fm(ig,k+1)+detr(ig,k))
    951950               else
    952951                  qa(ig,k)=q(ig,k)
    953                endif
     952               ENDIF
    954953               
    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
     954!               IF (qa(ig,k).lt.0.) print *, 'WARNING: qa is negative!'
     955!               IF (q(ig,k).lt.0.) print *, 'WARNING: q is negative!'
     956            ENDDO
     957         ENDDO
    959958         
    960959!------------------------------------------------------------------------------
    … …  
    965964            do ig=1,ngrid
    966965               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
     966               IF (wqd(ig,k).lt.0.) print*,'WARNING: wqd is negative!'
     967            ENDDO
     968         ENDDO
    970969         
    971970         do ig=1,ngrid
    972971            wqd(ig,1) = 0.
    973972            wqd(ig,nlay+1) = 0.
    974          enddo
     973         ENDDO
    975974         
    976975!------------------------------------------------------------------------------
    … …  
    983982               &       * (detr(ig,k) * qa(ig,k) - entr(ig,k) * q(ig,k)        &
    984983               &       - wqd(ig,k) + wqd(ig,k+1))
    985             enddo
    986          enddo
    987       endif
     984            ENDDO
     985         ENDDO
     986      ENDIF
    988987     
    989988      therm_tke_max(:,:) = q(:,:)
Note: See TracChangeset for help on using the changeset viewer.