! ! $Id: thermcell_main.F90 2193 2015-02-04 06:48:02Z musat $ ! SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep & & ,pplay,pplev,pphi,debut & & ,pu,pv,pt,po & & ,pduadj,pdvadj,pdtadj,pdoadj & & ,fm0,entr0,detr0,zqta,zqla,lmax & & ,ratqscth,ratqsdiff,zqsatth & & ,Ale_bl,Alp_bl,lalim_conv,wght_th & & ,zmax0, f0,zw2,fraca,ztv & & ,zpspsk,ztla,zthl & !!! nrlmd le 10/04/2012 & ,pbl_tke,pctsrf,omega,airephy & & ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 & & ,n2,s2,ale_bl_stat & & ,therm_tke_max,env_tke_max & & ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke & & ,alp_bl_conv,alp_bl_stat & !!! fin nrlmd le 10/04/2012 & ,ztva ) USE dimphy USE ioipsl USE comgeomphy , ONLY:rlond,rlatd USE indice_sol_mod IMPLICIT NONE !======================================================================= ! Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu ! Version du 09.02.07 ! Calcul du transport vertical dans la couche limite en presence ! de "thermiques" explicitement representes avec processus nuageux ! ! Reecriture a partir d'un listing papier a Habas, le 14/02/00 ! ! le thermique est suppose homogene et dissipe par melange avec ! son environnement. la longueur l_mix controle l'efficacite du ! melange ! ! Le calcul du transport des differentes especes se fait en prenant ! en compte: ! 1. un flux de masse montant ! 2. un flux de masse descendant ! 3. un entrainement ! 4. un detrainement ! ! Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr) ! Introduction of an implicit computation of vertical advection in ! the environment of thermal plumes in thermcell_dq ! impl = 0 : explicit, 1 : implicit, -1 : old version ! controled by iflag_thermals = ! 15, 16 run with impl=-1 : numerical convergence with NPv3 ! 17, 18 run with impl=1 : more stable ! 15 and 17 correspond to the activation of the stratocumulus "bidouille" ! !======================================================================= !----------------------------------------------------------------------- ! declarations: ! ------------- #include "dimensions.h" #include "YOMCST.h" #include "YOETHF.h" #include "FCTTRE.h" #include "iniprint.h" #include "thermcell.h" ! arguments: ! ---------- !IM 140508 INTEGER itap INTEGER ngrid,nlay real ptimestep REAL pt(ngrid,nlay),pdtadj(ngrid,nlay) REAL pu(ngrid,nlay),pduadj(ngrid,nlay) REAL pv(ngrid,nlay),pdvadj(ngrid,nlay) REAL po(ngrid,nlay),pdoadj(ngrid,nlay) REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1) real pphi(ngrid,nlay) LOGICAL debut ! local: ! ------ integer icount integer, save :: dvdq=1,dqimpl=-1 !$OMP THREADPRIVATE(dvdq,dqimpl) data icount/0/ save icount !$OMP THREADPRIVATE(icount) integer,save :: igout=1 !$OMP THREADPRIVATE(igout) integer,save :: lunout1=6 !$OMP THREADPRIVATE(lunout1) integer,save :: lev_out=10 !$OMP THREADPRIVATE(lev_out) REAL susqr2pi, Reuler INTEGER ig,k,l,ll,ierr real zsortie1d(klon) INTEGER lmax(klon),lmin(klon),lalim(klon) INTEGER lmix(klon) INTEGER lmix_bis(klon) real linter(klon) real zmix(klon) real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev) ! real fraca(klon,klev) real zmax_sec(klon) !on garde le zmax du pas de temps precedent real zmax0(klon) !FH/IM save zmax0 real lambda real zlev(klon,klev+1),zlay(klon,klev) real deltaz(klon,klev) REAL zh(klon,klev) real zthl(klon,klev),zdthladj(klon,klev) REAL ztv(klon,klev) real zu(klon,klev),zv(klon,klev),zo(klon,klev) real zl(klon,klev) real zsortie(klon,klev) real zva(klon,klev) real zua(klon,klev) real zoa(klon,klev) real zta(klon,klev) real zha(klon,klev) real fraca(klon,klev+1) real zf,zf2 real thetath2(klon,klev),wth2(klon,klev),wth3(klon,klev) real q2(klon,klev) ! FH probleme de dimensionnement avec l'allocation dynamique ! common/comtherm/thetath2,wth2 real wq(klon,klev) real wthl(klon,klev) real wthv(klon,klev) real ratqscth(klon,klev) real var real vardiff real ratqsdiff(klon,klev) logical sorties real rho(klon,klev),rhobarz(klon,klev),masse(klon,klev) real zpspsk(klon,klev) real wmax(klon) real wmax_tmp(klon) real wmax_sec(klon) real fm0(klon,klev+1),entr0(klon,klev),detr0(klon,klev) real fm(klon,klev+1),entr(klon,klev),detr(klon,klev) real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev) !niveau de condensation integer nivcon(klon) real zcon(klon) REAL CHI real zcon2(klon) real pcon(klon) real zqsat(klon,klev) real zqsatth(klon,klev) real f_star(klon,klev+1),entr_star(klon,klev) real detr_star(klon,klev) real alim_star_tot(klon) real alim_star(klon,klev) real alim_star_clos(klon,klev) real f(klon), f0(klon) !FH/IM save f0 real zlevinter(klon) real seuil real csc(klon,klev) !!! nrlmd le 10/04/2012 !------Entrées real pbl_tke(klon,klev+1,nbsrf) real pctsrf(klon,nbsrf) real omega(klon,klev) real airephy(klon) !------Sorties real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon) real therm_tke_max0(klon),env_tke_max0(klon) real n2(klon),s2(klon) real ale_bl_stat(klon) real therm_tke_max(klon,klev),env_tke_max(klon,klev) real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon) !------Local integer nsrf real rhobarz0(klon) ! Densité au LCL logical ok_lcl(klon) ! Existence du LCL des thermiques integer klcl(klon) ! Niveau du LCL real interp(klon) ! Coef d'interpolation pour le LCL !--Triggering real Su ! Surface unité: celle d'un updraft élémentaire parameter(Su=4e4) real hcoef ! Coefficient directeur pour le calcul de s2 parameter(hcoef=1) real hmincoef ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2 parameter(hmincoef=0.3) real eps1 ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd) parameter(eps1=0.3) real hmin(ngrid) ! Ordonnée à l'origine pour le calcul de s2 real zmax_moy(ngrid) ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl) real zmax_moy_coef parameter(zmax_moy_coef=0.33) real depth(klon) ! Epaisseur moyenne du cumulus real w_max(klon) ! Vitesse max statistique real s_max(klon) !--Closure real pbl_tke_max(klon,klev) ! Profil de TKE moyenne real pbl_tke_max0(klon) ! TKE moyenne au LCL real w_ls(klon,klev) ! Vitesse verticale grande échelle (m/s) real coef_m ! On considère un rendement pour alp_bl_fluct_m parameter(coef_m=1.) real coef_tke ! On considère un rendement pour alp_bl_fluct_tke parameter(coef_tke=1.) !!! fin nrlmd le 10/04/2012 ! !nouvelles variables pour la convection real Ale_bl(klon) real Alp_bl(klon) real alp_int(klon),dp_int(klon),zdp real ale_int(klon) integer n_int(klon) real fm_tot(klon) real wght_th(klon,klev) integer lalim_conv(klon) !v1d logical therm !v1d save therm character*2 str2 character*10 str10 character (len=20) :: modname='thermcell_main' character (len=80) :: abort_message EXTERNAL SCOPY ! !----------------------------------------------------------------------- ! initialisation: ! --------------- ! seuil=0.25 if (debut) then if (iflag_thermals==15.or.iflag_thermals==16) then dvdq=0 dqimpl=-1 else dvdq=1 dqimpl=1 endif fm0=0. entr0=0. detr0=0. endif fm=0. ; entr=0. ; detr=0. icount=icount+1 !IM 090508 beg !print*,'=====================================================================' !print*,'=====================================================================' !print*,' PAS ',icount,' PAS ',icount,' PAS ',icount,' PAS ',icount !print*,'=====================================================================' !print*,'=====================================================================' !IM 090508 end if (prt_level.ge.1) print*,'thermcell_main V4' sorties=.true. IF(ngrid.NE.klon) THEN PRINT* PRINT*,'STOP dans convadj' PRINT*,'ngrid =',ngrid PRINT*,'klon =',klon ENDIF ! ! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)' do ig=1,klon f0(ig)=max(f0(ig),1.e-2) zmax0(ig)=max(zmax0(ig),40.) !IMmarche pas ?! if (f0(ig)<1.e-2) f0(ig)=1.e-2 enddo if (prt_level.ge.20) then do ig=1,ngrid print*,'th_main ig f0',ig,f0(ig) enddo endif !----------------------------------------------------------------------- ! Calcul de T,q,ql a partir de Tl et qT dans l environnement ! -------------------------------------------------------------------- ! CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay, & & pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out) if (prt_level.ge.1) print*,'thermcell_main apres thermcell_env' !------------------------------------------------------------------------ ! -------------------- ! ! ! + + + + + + + + + + + ! ! ! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz ! wh,wt,wo ... ! ! + + + + + + + + + + + zh,zu,zv,zo,rho ! ! ! -------------------- zlev(1) ! \\\\\\\\\\\\\\\\\\\\ ! ! !----------------------------------------------------------------------- ! Calcul des altitudes des couches !----------------------------------------------------------------------- do l=2,nlay zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG enddo zlev(:,1)=0. zlev(:,nlay+1)=(2.*pphi(:,klev)-pphi(:,klev-1))/RG do l=1,nlay zlay(:,l)=pphi(:,l)/RG enddo !calcul de l epaisseur des couches do l=1,nlay deltaz(:,l)=zlev(:,l+1)-zlev(:,l) enddo ! print*,'2 OK convect8' !----------------------------------------------------------------------- ! Calcul des densites !----------------------------------------------------------------------- rho(:,:)=pplay(:,:)/(zpspsk(:,:)*RD*ztv(:,:)) if (prt_level.ge.10)write(lunout,*) & & 'WARNING thermcell_main rhobarz(:,1)=rho(:,1)' rhobarz(:,1)=rho(:,1) do l=2,nlay rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1)) enddo !calcul de la masse do l=1,nlay masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG enddo if (prt_level.ge.1) print*,'thermcell_main apres initialisation' !------------------------------------------------------------------ ! ! /|\ ! -------- | F_k+1 ------- ! ----> D_k ! /|\ <---- E_k , A_k ! -------- | F_k --------- ! ----> D_k-1 ! <---- E_k-1 , A_k-1 ! ! ! ! ! ! --------------------------- ! ! ----- F_lmax+1=0 ---------- \ ! lmax (zmax) | ! --------------------------- | ! | ! --------------------------- | ! | ! --------------------------- | ! | ! --------------------------- | ! | ! --------------------------- | ! | E ! --------------------------- | D ! | ! --------------------------- | ! | ! --------------------------- \ | ! lalim | | ! --------------------------- | | ! | | ! --------------------------- | | ! | A | ! --------------------------- | | ! | | ! --------------------------- | | ! lmin (=1 pour le moment) | | ! ----- F_lmin=0 ------------ / / ! ! --------------------------- ! ////////////////////////// ! ! !============================================================================= ! Calculs initiaux ne faisant pas intervenir les changements de phase !============================================================================= !------------------------------------------------------------------ ! 1. alim_star est le profil vertical de l'alimentation a la base du ! panache thermique, calcule a partir de la flotabilite de l'air sec ! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star !------------------------------------------------------------------ ! entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0. lmin=1 !----------------------------------------------------------------------------- ! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un ! panache sec conservatif (e=d=0) alimente selon alim_star ! Il s'agit d'un calcul de type CAPE ! zmax_sec est utilise pour determiner la geometrie du thermique. !------------------------------------------------------------------------------ !--------------------------------------------------------------------------------- !calcul du melange et des variables dans le thermique !-------------------------------------------------------------------------------- ! if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out !IM 140508 CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz, & ! Gestion temporaire de plusieurs appels à thermcell_plume au travers ! de la variable iflag_thermals ! print*,'THERM thermcell_main iflag_thermals_ed=',iflag_thermals_ed if (iflag_thermals_ed<=9) then ! print*,'THERM NOUVELLE/NOUVELLE Arnaud' CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,& & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & & ,lev_out,lunout1,igout) elseif (iflag_thermals_ed>9) then ! print*,'THERM RIO et al 2010, version d Arnaud' CALL thermcellV1_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,& & zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot, & & lalim,f0,detr_star,entr_star,f_star,csc,ztva, & & ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter & & ,lev_out,lunout1,igout) endif if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ') call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix ') if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume' if (prt_level.ge.10) then write(lunout1,*) 'Dans thermcell_main 2' write(lunout1,*) 'lmin ',lmin(igout) write(lunout1,*) 'lalim ',lalim(igout) write(lunout1,*) ' ig l alim_star entr_star detr_star f_star ' write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) & & ,f_star(igout,l+1),l=1,nint(linter(igout))+5) endif !------------------------------------------------------------------------------- ! Calcul des caracteristiques du thermique:zmax,zmix,wmax !------------------------------------------------------------------------------- ! CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2, & & zlev,lmax,zmax,zmax0,zmix,wmax,lev_out) ! Attention, w2 est transforme en sa racine carree dans cette routine ! Le probleme vient du fait que linter et lmix sont souvent égaux à 1. wmax_tmp=0. do l=1,nlay wmax_tmp(:)=max(wmax_tmp(:),zw2(:,l)) enddo ! print*,"ZMAX ",lalim,lmin,linter,lmix,lmax,zmax,zmax0,zmix,wmax call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ') call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ') call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ') call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax ') if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height' !------------------------------------------------------------------------------- ! Fermeture,determination de f !------------------------------------------------------------------------------- ! ! !! write(lunout,*)'THERM NOUVEAU XXXXX' CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star, & & lalim,lmin,zmax_sec,wmax_sec,lev_out) call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ') call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lalim ') if (prt_level.ge.1) print*,'thermcell_main apres thermcell_dry' if (prt_level.ge.10) then write(lunout1,*) 'Dans thermcell_main 1b' write(lunout1,*) 'lmin ',lmin(igout) write(lunout1,*) 'lalim ',lalim(igout) write(lunout1,*) ' ig l alim_star entr_star detr_star f_star ' write(lunout1,'(i6,i4,e15.5)') (igout,l,alim_star(igout,l) & & ,l=1,lalim(igout)+4) endif ! Choix de la fonction d'alimentation utilisee pour la fermeture. ! Apparemment sans importance alim_star_clos(:,:)=alim_star(:,:) alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:) ! !CR Appel de la fermeture seche if (iflag_thermals_closure.eq.1) then CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, & & zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Appel avec les zmax et wmax tenant compte de la condensation ! Semble moins bien marcher else if (iflag_thermals_closure.eq.2) then CALL thermcell_closure(ngrid,nlay,r_aspect_thermals,ptimestep,rho, & & zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out) endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure' if (tau_thermals>1.) then lambda=exp(-ptimestep/tau_thermals) f0=(1.-lambda)*f+lambda*f0 else f0=f endif ! Test valable seulement en 1D mais pas genant if (.not. (f0(1).ge.0.) ) then abort_message = '.not. (f0(1).ge.0.)' CALL abort_gcm (modname,abort_message,1) endif !------------------------------------------------------------------------------- !deduction des flux !------------------------------------------------------------------------------- CALL thermcell_flux2(ngrid,nlay,ptimestep,masse, & & lalim,lmax,alim_star, & & entr_star,detr_star,f,rhobarz,zlev,zw2,fm,entr, & & detr,zqla,lev_out,lunout1,igout) !IM 060508 & detr,zqla,zmax,lev_out,lunout,igout) if (prt_level.ge.1) print*,'thermcell_main apres thermcell_flux' call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ') call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax ') !------------------------------------------------------------------ ! On ne prend pas directement les profils issus des calculs precedents ! mais on s'autorise genereusement une relaxation vers ceci avec ! une constante de temps tau_thermals (typiquement 1800s). !------------------------------------------------------------------ if (tau_thermals>1.) then lambda=exp(-ptimestep/tau_thermals) fm0=(1.-lambda)*fm+lambda*fm0 entr0=(1.-lambda)*entr+lambda*entr0 detr0=(1.-lambda)*detr+lambda*detr0 else fm0=fm entr0=entr detr0=detr endif !c------------------------------------------------------------------ ! calcul du transport vertical !------------------------------------------------------------------ call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, & & zthl,zdthladj,zta,lev_out) call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, & & po,pdoadj,zoa,lev_out) !------------------------------------------------------------------ ! Calcul de la fraction de l'ascendance !------------------------------------------------------------------ do ig=1,klon fraca(ig,1)=0. fraca(ig,nlay+1)=0. enddo do l=2,nlay do ig=1,klon if (zw2(ig,l).gt.1.e-10) then fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l)) else fraca(ig,l)=0. endif enddo enddo !------------------------------------------------------------------ ! calcul du transport vertical du moment horizontal !------------------------------------------------------------------ !IM 090508 if (dvdq == 0 ) then ! Calcul du transport de V tenant compte d'echange par gradient ! de pression horizontal avec l'environnement call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse & ! & ,fraca*dvdq,zmax & & ,fraca,zmax & & ,zu,zv,pduadj,pdvadj,zua,zva,lev_out) else ! calcul purement conservatif pour le transport de V call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse & & ,zu,pduadj,zua,lev_out) call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse & & ,zv,pdvadj,zva,lev_out) endif ! print*,'13 OK convect8' do l=1,nlay do ig=1,ngrid pdtadj(ig,l)=zdthladj(ig,l)*zpspsk(ig,l) enddo enddo if (prt_level.ge.1) print*,'14 OK convect8' !------------------------------------------------------------------ ! Calculs de diagnostiques pour les sorties !------------------------------------------------------------------ !calcul de fraca pour les sorties if (sorties) then if (prt_level.ge.1) print*,'14a OK convect8' ! calcul du niveau de condensation ! initialisation do ig=1,ngrid nivcon(ig)=0 zcon(ig)=0. enddo !nouveau calcul do ig=1,ngrid CHI=zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1)) pcon(ig)=pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI enddo !IM do k=1,nlay do k=1,nlay-1 do ig=1,ngrid if ((pcon(ig).le.pplay(ig,k)) & & .and.(pcon(ig).gt.pplay(ig,k+1))) then zcon2(ig)=zlay(ig,k)-(pcon(ig)-pplay(ig,k))/(RG*rho(ig,k))/100. endif enddo enddo !IM ierr=0 do ig=1,ngrid if (pcon(ig).le.pplay(ig,nlay)) then zcon2(ig)=zlay(ig,nlay)-(pcon(ig)-pplay(ig,nlay))/(RG*rho(ig,nlay))/100. ierr=1 endif enddo if (ierr==1) then abort_message = 'thermcellV0_main: les thermiques vont trop haut ' CALL abort_gcm (modname,abort_message,1) endif if (prt_level.ge.1) print*,'14b OK convect8' do k=nlay,1,-1 do ig=1,ngrid if (zqla(ig,k).gt.1e-10) then nivcon(ig)=k zcon(ig)=zlev(ig,k) endif enddo enddo if (prt_level.ge.1) print*,'14c OK convect8' !calcul des moments !initialisation do l=1,nlay do ig=1,ngrid q2(ig,l)=0. wth2(ig,l)=0. wth3(ig,l)=0. ratqscth(ig,l)=0. ratqsdiff(ig,l)=0. enddo enddo if (prt_level.ge.1) print*,'14d OK convect8' if (prt_level.ge.10)write(lunout,*) & & 'WARNING thermcell_main wth2=0. si zw2 > 1.e-10' do l=1,nlay do ig=1,ngrid zf=fraca(ig,l) zf2=zf/(1.-zf) ! thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2 if(zw2(ig,l).gt.1.e-10) then wth2(ig,l)=zf2*(zw2(ig,l))**2 else wth2(ig,l)=0. endif wth3(ig,l)=zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l)) & & *zw2(ig,l)*zw2(ig,l)*zw2(ig,l) q2(ig,l)=zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2 !test: on calcul q2/po=ratqsc ratqscth(ig,l)=sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.)) enddo enddo !calcul des flux: q, thetal et thetav do l=1,nlay do ig=1,ngrid wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.) wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l)) wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l)) enddo enddo ! !!! nrlmd le 10/04/2012 !------------Test sur le LCL des thermiques do ig=1,ngrid ok_lcl(ig)=.false. if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true. enddo !------------Localisation des niveaux entourant le LCL et du coef d'interpolation do l=1,nlay-1 do ig=1,ngrid if (ok_lcl(ig)) then !ATTENTION,zw2 calcule en pplev ! if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then ! klcl(ig)=l ! interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig))) ! endif if ((pplev(ig,l) .ge. pcon(ig)) .and. (pplev(ig,l+1) .le. pcon(ig))) then klcl(ig)=l interp(ig)=(pcon(ig)-pplev(ig,klcl(ig)))/(pplev(ig,klcl(ig)+1)-pplev(ig,klcl(ig))) endif endif enddo enddo !------------Hauteur des thermiques !!jyg le 27/04/2012 !! do ig =1,ngrid !! rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) & !! & -rhobarz(ig,klcl(ig)))*interp(ig) !! zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG) !! zmax(ig)=pphi(ig,lmax(ig))/rg !! if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax !! enddo do ig =1,ngrid !CR:REHABILITATION ZMAX CONTINU ! zmax(ig)=pphi(ig,lmax(ig))/rg if (ok_lcl(ig)) then rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) & & -rhobarz(ig,klcl(ig)))*interp(ig) zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG) zlcl(ig)=min(zlcl(ig),zmax(ig)) ! Si zlcl > zmax alors on pose zlcl = zmax else rhobarz0(ig)=0. zlcl(ig)=zmax(ig) endif enddo !!jyg fin !------------Calcul des propriétés du thermique au LCL IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN !-----Initialisation de la TKE moyenne do l=1,nlay do ig=1,ngrid pbl_tke_max(ig,l)=0. enddo enddo !-----Calcul de la TKE moyenne do nsrf=1,nbsrf do l=1,nlay do ig=1,ngrid pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l) enddo enddo enddo !-----Initialisations des TKE dans et hors des thermiques do l=1,nlay do ig=1,ngrid therm_tke_max(ig,l)=pbl_tke_max(ig,l) env_tke_max(ig,l)=pbl_tke_max(ig,l) enddo enddo !-----Calcul de la TKE transportée par les thermiques : therm_tke_max call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, & & rg,pplev,therm_tke_max) ! print *,' thermcell_tke_transport -> ' !!jyg !-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls do l=1,nlay do ig=1,ngrid pbl_tke_max(ig,l)=fraca(ig,l)*therm_tke_max(ig,l)+(1.-fraca(ig,l))*env_tke_max(ig,l) ! Recalcul de TKE moyenne aprés transport de TKE_TH env_tke_max(ig,l)=(pbl_tke_max(ig,l)-fraca(ig,l)*therm_tke_max(ig,l))/(1.-fraca(ig,l)) ! Recalcul de TKE dans l'environnement aprés transport de TKE_TH w_ls(ig,l)=-1.*omega(ig,l)/(RG*rhobarz(ig,l)) ! Vitesse verticale de grande échelle enddo enddo ! print *,' apres w_ls = ' !!jyg do ig=1,ngrid if (ok_lcl(ig)) then fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) & & -fraca(ig,klcl(ig)))*interp(ig) w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) & & -zw2(ig,klcl(ig)))*interp(ig) w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) & & -w_ls(ig,klcl(ig)))*interp(ig) therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) & & +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig) env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) & & -env_tke_max(ig,klcl(ig)))*interp(ig) pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) & & -pbl_tke_max(ig,klcl(ig)))*interp(ig) if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20. if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20. if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20. else fraca0(ig)=0. w0(ig)=0. !!jyg le 27/04/2012 !! zlcl(ig)=0. !! endif enddo ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ! print *,'ENDIF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) ' !!jyg !------------Triggering------------------ IF (iflag_trig_bl.ge.1) THEN !-----Initialisations depth(:)=0. n2(:)=0. s2(:)=100. ! some low value, arbitrary s_max(:)=0. !-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max) do ig=1,ngrid zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig)) depth(ig)=zmax_moy(ig)-zlcl(ig) hmin(ig)=hmincoef*zlcl(ig) if (depth(ig).ge.10.) then s2(ig)=(hcoef*depth(ig)+hmin(ig))**2 n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig) !! !!jyg le 27/04/2012 !! s_max(ig)=s2(ig)*log(n2(ig)) !! if (n2(ig) .lt. 1) s_max(ig)=0. s_max(ig)=s2(ig)*log(max(n2(ig),1.)) !!fin jyg else n2(ig)=0. s_max(ig)=0. endif enddo ! print *,'avant Calcul de Wmax ' !!jyg !-----Calcul de Wmax et ALE_BL_STAT associée !!jyg le 30/04/2012 !! do ig=1,ngrid !! if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.1.) ) then !! w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/su)-log(2.*3.14)-log(2.*log(s_max(ig)/su)-log(2.*3.14)))) !! ale_bl_stat(ig)=0.5*w_max(ig)**2 !! else !! w_max(ig)=0. !! ale_bl_stat(ig)=0. !! endif !! enddo susqr2pi=su*sqrt(2.*Rpi) Reuler=exp(1.) do ig=1,ngrid if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*Reuler) ) then w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi)))) ale_bl_stat(ig)=0.5*w_max(ig)**2 else w_max(ig)=0. ale_bl_stat(ig)=0. endif enddo ENDIF ! iflag_trig_bl ! print *,'ENDIF iflag_trig_bl' !!jyg !------------Closure------------------ IF (iflag_clos_bl.ge.2) THEN !-----Calcul de ALP_BL_STAT do ig=1,ngrid alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2) alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* & & (w0(ig)**2) alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) & & +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig) if (iflag_clos_bl.ge.2) then alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* & & (w0(ig)**2) else alp_bl_conv(ig)=0. endif alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig) enddo !-----Sécurité ALP infinie do ig=1,ngrid if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2. enddo ENDIF ! (iflag_clos_bl.ge.2) !!! fin nrlmd le 10/04/2012 if (prt_level.ge.10) then ig=igout do l=1,nlay print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l) print*,'14g OK convect8 ig,l,po',ig,l,po(ig,l) enddo endif ! print*,'avant calcul ale et alp' !calcul de ALE et ALP pour la convection Alp_bl(:)=0. Ale_bl(:)=0. ! print*,'ALE,ALP ,l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig)' do l=1,nlay do ig=1,ngrid Alp_bl(ig)=max(Alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) ) Ale_bl(ig)=max(Ale_bl(ig),0.5*zw2(ig,l)**2) ! print*,'ALE,ALP',l,zw2(ig,l),Ale_bl(ig),Alp_bl(ig) enddo enddo ! Ale sec (max de wmax/2 sous la zone d'inhibition) dans ! le cas iflag_trig_bl=3 IF (iflag_trig_bl==3) Ale_bl(:)=0.5*wmax_sec(:)**2 !test:calcul de la ponderation des couches pour KE !initialisations fm_tot(:)=0. wght_th(:,:)=1. lalim_conv(:)=lalim(:) do k=1,klev do ig=1,ngrid if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k) enddo enddo ! assez bizarre car, si on est dans la couche d'alim et que alim_star et ! plus petit que 1.e-10, on prend wght_th=1. do k=1,klev do ig=1,ngrid if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then wght_th(ig,k)=alim_star(ig,k) endif enddo enddo ! print*,'apres wght_th' !test pour prolonger la convection do ig=1,ngrid !v1d if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then if ((alim_star(ig,1).lt.1.e-10)) then lalim_conv(ig)=1 wght_th(ig,1)=1. ! print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1) endif enddo !------------------------------------------------------------------------ ! Modif CR/FH 20110310 : Alp integree sur la verticale. ! Integrale verticale de ALP. ! wth3 etant aux niveaux inter-couches, on utilise d play comme masse des ! couches !------------------------------------------------------------------------ alp_int(:)=0. dp_int(:)=0. do l=2,nlay do ig=1,ngrid if(l.LE.lmax(ig)) THEN zdp=pplay(ig,l-1)-pplay(ig,l) alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)*zdp dp_int(ig)=dp_int(ig)+zdp endif enddo enddo if (iflag_coupl>=3 .and. iflag_coupl<=5) then do ig=1,ngrid !valeur integree de alp_bl * 0.5: if (dp_int(ig)>0.) then Alp_bl(ig)=alp_int(ig)/dp_int(ig) endif enddo! endif ! Facteur multiplicatif sur Alp_bl Alp_bl(:)=alp_bl_k*Alp_bl(:) !------------------------------------------------------------------------ !calcul du ratqscdiff if (prt_level.ge.1) print*,'14e OK convect8' var=0. vardiff=0. ratqsdiff(:,:)=0. do l=1,klev do ig=1,ngrid if (l<=lalim(ig)) then var=var+alim_star(ig,l)*zqta(ig,l)*1000. endif enddo enddo if (prt_level.ge.1) print*,'14f OK convect8' do l=1,klev do ig=1,ngrid if (l<=lalim(ig)) then zf=fraca(ig,l) zf2=zf/(1.-zf) vardiff=vardiff+alim_star(ig,l)*(zqta(ig,l)*1000.-var)**2 endif enddo enddo if (prt_level.ge.1) print*,'14g OK convect8' do l=1,nlay do ig=1,ngrid ratqsdiff(ig,l)=sqrt(vardiff)/(po(ig,l)*1000.) ! write(11,*)'ratqsdiff=',ratqsdiff(ig,l) enddo enddo !-------------------------------------------------------------------- ! !ecriture des fichiers sortie ! print*,'15 OK convect8 CCCCCCCCCCCCCCCCCCc' endif if (prt_level.ge.1) print*,'thermcell_main FIN OK' return end !----------------------------------------------------------------------------- subroutine test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,ztva,zqla,f_star,zw2,comment) IMPLICIT NONE #include "iniprint.h" integer i, k, klon,klev real pplev(klon,klev+1),pplay(klon,klev) real ztv(klon,klev) real po(klon,klev) real ztva(klon,klev) real zqla(klon,klev) real f_star(klon,klev) real zw2(klon,klev) integer long(klon) real seuil character*21 comment if (prt_level.ge.1) THEN print*,'WARNING !!! TEST ',comment endif return ! test sur la hauteur des thermiques ... do i=1,klon !IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then if (prt_level.ge.10) then print*,'WARNING ',comment,' au point ',i,' K= ',long(i) print*,' K P(MB) THV(K) Qenv(g/kg)THVA QLA(g/kg) F* W2' do k=1,klev 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) enddo endif enddo return end !!! nrlmd le 10/04/2012 Transport de la TKE par le thermique moyen pour la fermeture en ALP ! On transporte pbl_tke pour donner therm_tke ! Copie conforme de la subroutine DTKE dans physiq.F écrite par Frederic Hourdin subroutine thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0, & & rg,pplev,therm_tke_max) implicit none #include "iniprint.h" !======================================================================= ! ! Calcul du transport verticale dans la couche limite en presence ! de "thermiques" explicitement representes ! calcul du dq/dt une fois qu'on connait les ascendances ! !======================================================================= integer ngrid,nlay,nsrf real ptimestep real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1) real entr0(ngrid,nlay),rg real therm_tke_max(ngrid,nlay) real detr0(ngrid,nlay) real masse(ngrid,nlay),fm(ngrid,nlay+1) real entr(ngrid,nlay) real q(ngrid,nlay) integer lev_out ! niveau pour les print real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1) real zzm integer ig,k integer isrf lev_out=0 if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0' ! calcul du detrainement do k=1,nlay detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k) masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG enddo ! Decalage vertical des entrainements et detrainements. masse(:,1)=0.5*masse0(:,1) entr(:,1)=0.5*entr0(:,1) detr(:,1)=0.5*detr0(:,1) fm(:,1)=0. do k=1,nlay-1 masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1)) entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1)) detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1)) fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k) enddo fm(:,nlay+1)=0. !!! nrlmd le 16/09/2010 ! calcul de la valeur dans les ascendances ! do ig=1,ngrid ! qa(ig,1)=q(ig,1) ! enddo !!! !do isrf=1,nsrf ! q(:,:)=therm_tke(:,:,isrf) q(:,:)=therm_tke_max(:,:) !!! nrlmd le 16/09/2010 do ig=1,ngrid qa(ig,1)=q(ig,1) enddo !!! if (1==1) then do k=2,nlay do ig=1,ngrid if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt. & & 1.e-5*masse(ig,k)) then qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k)) & & /(fm(ig,k+1)+detr(ig,k)) else qa(ig,k)=q(ig,k) endif if (qa(ig,k).lt.0.) then ! print*,'qa<0!!!' endif if (q(ig,k).lt.0.) then ! print*,'q<0!!!' endif enddo enddo ! Calcul du flux subsident do k=2,nlay do ig=1,ngrid wqd(ig,k)=fm(ig,k)*q(ig,k) if (wqd(ig,k).lt.0.) then ! print*,'wqd<0!!!' endif enddo enddo do ig=1,ngrid wqd(ig,1)=0. wqd(ig,nlay+1)=0. enddo ! Calcul des tendances do k=1,nlay do ig=1,ngrid q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k) & & -wqd(ig,k)+wqd(ig,k+1)) & & *ptimestep/masse(ig,k) enddo enddo endif therm_tke_max(:,:)=q(:,:) return !!! fin nrlmd le 10/04/2012 end