Index: LMDZ6/trunk/libf/phylmdiso/calltherm.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/calltherm.F90 (revision 4090) +++ (revision ) @@ -1,539 +1,0 @@ -! -! $Id$ -! - subroutine calltherm(dtime & - & ,pplay,paprs,pphi,weak_inversion & - & ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & - & ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs & - & ,fm_therm,entr_therm,detr_therm,zqasc,clwcon0,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_th,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 - & ,zqla,ztva & -#ifdef ISO - & ,xt_seri,d_xt_ajs & -#ifdef DIAGISO - & ,q_the,xt_the & -#endif -#endif - & ) - - USE dimphy - USE indice_sol_mod - USE print_control_mod, ONLY: prt_level,lunout -#ifdef ISO - use infotrac_phy, ONLY: ntraciso -#ifdef ISOVERIF - USE isotopes_mod, ONLY: iso_eau,iso_HDO - USE isotopes_verif_mod, ONLY: iso_verif_aberrant_enc_vect2D, & - iso_verif_egalite_vect2D -#endif -#endif - - implicit none - include "clesphys.h" - include "thermcell_old.h" - - -!IM 140508 - INTEGER, SAVE :: itap -!$OMP THREADPRIVATE(itap) - REAL dtime - LOGICAL debut - LOGICAL logexpr0, logexpr2(klon,klev), logexpr1(klon) - REAL fact(klon) - INTEGER nbptspb - - REAL u_seri(klon,klev),v_seri(klon,klev) - REAL t_seri(klon,klev),q_seri(klon,klev),qmemoire(klon,klev) - REAL weak_inversion(klon) - REAL paprs(klon,klev+1) - REAL pplay(klon,klev) - REAL pphi(klon,klev) - real zlev(klon,klev+1) -!test: on sort lentr et a* pour alimenter KE - REAL wght_th(klon,klev) - INTEGER lalim_conv(klon) - REAL zw2(klon,klev+1),fraca(klon,klev+1) - -!FH Update Thermiques - REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev) - REAL d_u_ajs(klon,klev),d_v_ajs(klon,klev) - real fm_therm(klon,klev+1) - real entr_therm(klon,klev),detr_therm(klon,klev) - -!******************************************************** -! declarations - LOGICAL flag_bidouille_stratocu - real fmc_therm(klon,klev+1),zqasc(klon,klev) - real zqla(klon,klev) - real zqta(klon,klev) - real ztv(klon,klev),ztva(klon,klev) - real zpspsk(klon,klev) - real ztla(klon,klev) - real zthl(klon,klev) - real wmax_sec(klon) - real zmax_sec(klon) - real f_sec(klon) - real detrc_therm(klon,klev) -! FH WARNING : il semble que ces save ne servent a rien -! save fmc_therm, detrc_therm - real clwcon0(klon,klev) - real zqsat(klon,klev) - real zw_sec(klon,klev+1) - integer lmix_sec(klon) - integer lmax(klon) - real ratqscth(klon,klev) - real ratqsdiff(klon,klev) - real zqsatth(klon,klev) -!nouvelles variables pour la convection - real ale_bl(klon) - real alp_bl(klon) - real ale(klon) - real alp(klon) -!RC - !on garde le zmax du pas de temps precedent - real zmax0(klon), f0(klon) - -!!! nrlmd le 10/04/2012 - real pbl_tke(klon,klev+1,nbsrf) - real pctsrf(klon,nbsrf) - real omega(klon,klev) - real airephy(klon) - real zlcl_th(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) -!!! fin nrlmd le 10/04/2012 - -!******************************************************** - - real, dimension(klon) :: pcon - real, dimension(klon,klev) :: rhobarz,wth3 - integer,dimension(klon) :: lalim - real, dimension(klon,klev+1) :: fm - real, dimension(klon,klev) :: alim_star - real, dimension(klon) :: zmax - - - - -! variables locales - REAL d_t_the(klon,klev), d_q_the(klon,klev) - REAL d_u_the(klon,klev),d_v_the(klon,klev) -! - real zfm_therm(klon,klev+1),zdt - real zentr_therm(klon,klev),zdetr_therm(klon,klev) -! FH A VERIFIER : SAVE INUTILES -! save zentr_therm,zfm_therm - - character (len=20) :: modname='calltherm' - character (len=80) :: abort_message - - integer i,k,isplit - logical, save :: first=.true. - logical :: new_thermcell - -#ifdef ISO - REAL xt_seri(ntraciso,klon,klev),xtmemoire(ntraciso,klon,klev) - REAL d_xt_ajs(ntraciso,klon,klev) - real d_xt_the(ntraciso,klon,klev) -#ifdef DIAGISO - real q_the(klon,klev) - real xt_the(ntraciso,klon,klev) -#endif - real qprec(klon,klev) - integer ixt -#endif - - -!$OMP THREADPRIVATE(first) -!******************************************************** - if (first) then - itap=0 - first=.false. - endif - -! Incrementer le compteur de la physique - itap = itap + 1 - -! Modele du thermique -! =================== -! print*,'thermiques: WARNING on passe t au lieu de t_seri' - - -! On prend comme valeur initiale des thermiques la valeur du pas -! de temps precedent - zfm_therm(:,:)=fm_therm(:,:) - zdetr_therm(:,:)=detr_therm(:,:) - zentr_therm(:,:)=entr_therm(:,:) - -! On reinitialise les flux de masse a zero pour le cumul en -! cas de splitting - fm_therm(:,:)=0. - entr_therm(:,:)=0. - detr_therm(:,:)=0. - - ale_bl(:)=0. - alp_bl(:)=0. - if (prt_level.ge.10) then - print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion' - endif - -! tests sur les valeurs negatives de l'eau - logexpr0=prt_level.ge.10 - nbptspb=0 - do k=1,klev - do i=1,klon -! Attention teste abderr 19-03-09 -! logexpr2(i,k)=.not.q_seri(i,k).ge.0. - logexpr2(i,k)=.not.q_seri(i,k).ge.1.e-15 - if (logexpr2(i,k)) then -#ifdef ISO - qprec(i,k)=q_seri(i,k) -#endif - q_seri(i,k)=1.e-15 - nbptspb=nbptspb+1 -#ifdef ISO - do ixt=1,ntraciso - xt_seri(ixt,i,k)=1.e-15*(xt_seri(ixt,i,k)/qprec(i,k)) - ! xt_seri(ixt,i,k)=1.e-15*(Rdefault(index_iso(ixt))) - enddo -#endif - endif -! if (logexpr0) & -! & print*,'WARN eau<0 avant therm i=',i,' k=',k & -! & ,' dq,q',d_q_the(i,k),q_seri(i,k) - enddo - enddo - if(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb - - - new_thermcell=iflag_thermals>=15.and.iflag_thermals<=18 -#ifdef ISO - if (.not.new_thermcell) then - CALL abort_gcm('calltherm 234','isos pas prevus ici',1) - endif -#ifdef ISOVERIF - if (iso_eau.gt.0) then - call iso_verif_egalite_vect2D( & - & xt_seri,q_seri, & - & 'calltherm 174',ntraciso,klon,klev) - endif !if (iso_eau.gt.0) then -#endif -#endif - zdt=dtime/REAL(nsplit_thermals) - - - do isplit=1,nsplit_thermals - - if (iflag_thermals>=1000) then - CALL thermcell_2002(klon,klev,zdt,iflag_thermals & - & ,pplay,paprs,pphi & - & ,u_seri,v_seri,t_seri,q_seri & - & ,d_u_the,d_v_the,d_t_the,d_q_the & - & ,zfm_therm,zentr_therm,fraca,zw2 & - & ,r_aspect_thermals,30.,w2di_thermals & - & ,tau_thermals) - else if (iflag_thermals.eq.2) then - CALL thermcell_sec(klon,klev,zdt & - & ,pplay,paprs,pphi,zlev & - & ,u_seri,v_seri,t_seri,q_seri & - & ,d_u_the,d_v_the,d_t_the,d_q_the & - & ,zfm_therm,zentr_therm & - & ,r_aspect_thermals,30.,w2di_thermals & - & ,tau_thermals) - else if (iflag_thermals.eq.3) then - CALL thermcell(klon,klev,zdt & - & ,pplay,paprs,pphi & - & ,u_seri,v_seri,t_seri,q_seri & - & ,d_u_the,d_v_the,d_t_the,d_q_the & - & ,zfm_therm,zentr_therm & - & ,r_aspect_thermals,l_mix_thermals,w2di_thermals & - & ,tau_thermals) - else if (iflag_thermals.eq.10) then - CALL thermcell_eau(klon,klev,zdt & - & ,pplay,paprs,pphi & - & ,u_seri,v_seri,t_seri,q_seri & - & ,d_u_the,d_v_the,d_t_the,d_q_the & - & ,zfm_therm,zentr_therm & - & ,r_aspect_thermals,l_mix_thermals,w2di_thermals & - & ,tau_thermals) - else if (iflag_thermals.eq.11) then - abort_message = 'cas non prevu dans calltherm' - CALL abort_physic (modname,abort_message,1) - else if (iflag_thermals.eq.12) then - CALL calcul_sec(klon,klev,zdt & - & ,pplay,paprs,pphi,zlev & - & ,u_seri,v_seri,t_seri,q_seri & - & ,zmax_sec,wmax_sec,zw_sec,lmix_sec & - & ,r_aspect_thermals,l_mix_thermals,w2di_thermals & - & ,tau_thermals) - else if (iflag_thermals==13.or.iflag_thermals==14) then - abort_message = 'thermcellV0_main enleve svn>2084' - CALL abort_physic (modname,abort_message,1) - else if (new_thermcell) then - CALL thermcell_main(itap,klon,klev,zdt & - & ,pplay,paprs,pphi,debut & - & ,u_seri,v_seri,t_seri,q_seri & - & ,d_u_the,d_v_the,d_t_the,d_q_the & - & ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax & - & ,ratqscth,ratqsdiff,zqsatth & - & ,zmax0,f0,zw2,fraca,ztv,zpspsk & - & ,ztla,zthl,ztva & - & ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax & -#ifdef ISO - & ,xt_seri,d_xt_the & -#endif - & ) - - CALL thermcell_alp(klon,klev,zdt & ! in - & ,pplay,paprs & ! in - & ,zfm_therm,zentr_therm,lmax & ! in - & ,pbl_tke,pctsrf,omega,airephy & ! in - & ,zw2,fraca & ! in - & ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax & ! in - & ,ale,alp,lalim_conv,wght_th & ! out - & ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &! out - & ,n2,s2,ale_bl_stat & ! out - & ,therm_tke_max,env_tke_max & ! out - & ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke & ! out - & ,alp_bl_conv,alp_bl_stat & ! out - & ) - - if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK' - else - abort_message = 'Cas des thermiques non prevu' - CALL abort_physic (modname,abort_message,1) - endif - -! Attention : les noms sont contre intuitif. -! flag_bidouille_stratocu est .true. si on ne fait pas de bidouille. -! Il aurait mieux valu avoir un nobidouille_stratocu -! Et pour simplifier : -! nobidouille_stratocu=.not.(iflag_thermals==13.or.iflag_thermals=15) -! Ce serait bien de changer, mai en prenant le temps de vérifier que ca -! fait bien ce qu'on croit. - - flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18 - -! Calcul a posteriori du niveau max des thermiques pour les schémas qui -! ne la sortent pas. - if (iflag_thermals<=12.or.iflag_thermals>=1000) then - lmax(:)=1 - do k=1,klev-1 - zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1) - enddo - do k=1,klev-1 - do i=1,klon - if (zfm_therm(i,k+1)>0.) lmax(i)=k - enddo - enddo - endif - - fact(:)=0. - DO i=1,klon - logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5 - IF(logexpr1(i)) fact(i)=1./REAL(nsplit_thermals) - ENDDO - - DO k=1,klev -! transformation de la derivee en tendance - d_t_the(:,k)=d_t_the(:,k)*dtime*fact(:) - d_u_the(:,k)=d_u_the(:,k)*dtime*fact(:) - d_v_the(:,k)=d_v_the(:,k)*dtime*fact(:) - d_q_the(:,k)=d_q_the(:,k)*dtime*fact(:) - fm_therm(:,k)=fm_therm(:,k) & - & +zfm_therm(:,k)*fact(:) - entr_therm(:,k)=entr_therm(:,k) & - & +zentr_therm(:,k)*fact(:) - detr_therm(:,k)=detr_therm(:,k) & - & +zdetr_therm(:,k)*fact(:) -#ifdef ISO - do ixt=1,ntraciso - d_xt_the(ixt,:,k)=d_xt_the(ixt,:,k)*dtime*fact(:) - enddo -#endif - ENDDO - fm_therm(:,klev+1)=0. - - - -! accumulation de la tendance - d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:) - d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:) - d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:) - d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_the(:,:) -#ifdef ISO - d_xt_ajs(:,:,:)=d_xt_ajs(:,:,:)+d_xt_the(:,:,:) -#endif - -! incrementation des variables meteo - t_seri(:,:) = t_seri(:,:) + d_t_the(:,:) - u_seri(:,:) = u_seri(:,:) + d_u_the(:,:) - v_seri(:,:) = v_seri(:,:) + d_v_the(:,:) - qmemoire(:,:)=q_seri(:,:) - q_seri(:,:) = q_seri(:,:) + d_q_the(:,:) -#ifdef ISO - xtmemoire(:,:,:)=xt_seri(:,:,:) - xt_seri(:,:,:) = xt_seri(:,:,:) + d_xt_the(:,:,:) -#ifdef ISOVERIF -! write(*,*) 'calltherm 350 tmp: ajout d_xt_the' - if (iso_HDO.gt.0) then -! i=479 -! k=4 -! write(*,*) 'xt_seri(iso_hdo,i,k),q_seri(i,k)=', & -! & xt_seri(iso_hdo,i,k),q_seri(i,k) -! write(*,*) 'd_xt_the(iso_hdo,i,k),d_q_the(i,k)=', & -! & d_xt_the(iso_hdo,i,k),d_q_the(i,k) - call iso_verif_aberrant_enc_vect2D( & - & xt_seri,q_seri, & - & 'calltherm 353, apres ajout d_xt_the',ntraciso,klon,klev) - endif -#endif -#endif - if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK' - - DO i=1,klon - fm_therm(i,klev+1)=0. - ale_bl(i)=ale_bl(i)+ale(i)/REAL(nsplit_thermals) -! write(22,*)'ALE CALLTHERM',ale_bl(i),ale(i) - alp_bl(i)=alp_bl(i)+alp(i)/REAL(nsplit_thermals) -! write(23,*)'ALP CALLTHERM',alp_bl(i),alp(i) - if(prt_level.GE.10) print*,'calltherm i alp_bl alp ale_bl ale',i,alp_bl(i),alp(i),ale_bl(i),ale(i) - ENDDO - -!IM 060508 marche pas comme cela !!! enddo ! isplit - -! tests sur les valeurs negatives de l'eau - nbptspb=0 - DO k = 1, klev - DO i = 1, klon - logexpr2(i,k)=.not.q_seri(i,k).ge.0. - if (logexpr2(i,k)) then - q_seri(i,k)=1.e-15 - nbptspb=nbptspb+1 -#ifdef ISO - do ixt=1,ntraciso - xt_seri(ixt,i,k)=1.e-15*(xtmemoire(ixt,i,k)/qmemoire(i,k)) - enddo -#endif -! if (prt_level.ge.10) then -! print*,'WARN eau<0 apres therm i=',i,' k=',k & -! & ,' dq,q',d_q_the(i,k),q_seri(i,k), & -! & 'fm=',zfm_therm(i,k),'entr=',entr_therm(i,k) - endif - ENDDO - ENDDO -#ifdef ISO -#ifdef ISOVERIF - if (iso_HDO.gt.0) then - call iso_verif_aberrant_enc_vect2D( & - & xt_seri,q_seri, & - & 'calltherm 393, apres bidouille q<0',ntraciso,klon,klev) - endif -#endif -#endif - - IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb -! tests sur les valeurs de la temperature - nbptspb=0 - DO k = 1, klev - DO i = 1, klon - logexpr2(i,k)=t_seri(i,k).lt.50..or.t_seri(i,k).gt.370. - if (logexpr2(i,k)) nbptspb=nbptspb+1 -! if ((t_seri(i,k).lt.50.) .or. & -! & (t_seri(i,k).gt.370.)) then -! print*,'WARN temp apres therm i=',i,' k=',k & -! & ,' t_seri',t_seri(i,k) -! CALL abort -! endif - ENDDO - ENDDO - IF(nbptspb.GT.0) print*,'Number of points with q_seri(i,k)<=0 ',nbptspb - enddo ! isplit - -! -!*************************************************************** -! calcul du flux ascencant conservatif -! print*,'<<<>> - REAL, ALLOCATABLE, SAVE :: falb_dif(:,:,:), falb_dir(:,:,:) - REAL, ALLOCATABLE, SAVE :: chl_con(:) -!$OMP THREADPRIVATE(falb_dir,falb_dif,chl_con) -!albedo SB <<< - - - REAL, ALLOCATABLE, SAVE :: rain_fall(:), snow_fall(:) -!$OMP THREADPRIVATE( rain_fall, snow_fall) - REAL, ALLOCATABLE, SAVE :: solsw(:), solswfdiff(:), sollw(:) -!$OMP THREADPRIVATE(solsw, solswfdiff, sollw) - REAL, ALLOCATABLE, SAVE :: radsol(:) -!$OMP THREADPRIVATE(radsol) - REAL, ALLOCATABLE, SAVE :: swradcorr(:) -!$OMP THREADPRIVATE(swradcorr) -#ifdef ISO - REAL,ALLOCATABLE,SAVE :: xtsol(:,:),fxtevap(:,:,:) -!$OMP THREADPRIVATE(xtsol,fxtevap) - REAL, ALLOCATABLE, SAVE :: xtrain_fall(:,:), xtsnow_fall(:,:) -!$OMP THREADPRIVATE(xtrain_fall,xtsnow_fall) -#endif - -!clesphy0 param physiq -! -! Parametres de l'Orographie a l'Echelle Sous-Maille (OESM): -! - REAL, ALLOCATABLE, SAVE :: zmea(:), zstd(:), zsig(:), zgam(:) -!$OMP THREADPRIVATE(zmea, zstd, zsig, zgam) - REAL, ALLOCATABLE, SAVE :: zthe(:), zpic(:), zval(:) -!$OMP THREADPRIVATE(zthe, zpic, zval) -! REAL tabcntr0(100) - REAL, ALLOCATABLE, SAVE :: rugoro(:) -!$OMP THREADPRIVATE(rugoro) - REAL, ALLOCATABLE, SAVE :: t_ancien(:,:), q_ancien(:,:) -!$OMP THREADPRIVATE(t_ancien, q_ancien) - REAL, ALLOCATABLE, SAVE :: ql_ancien(:,:), qs_ancien(:,:) -!$OMP THREADPRIVATE(ql_ancien, qs_ancien) - REAL, ALLOCATABLE, SAVE :: prw_ancien(:), prlw_ancien(:), prsw_ancien(:) -!$OMP THREADPRIVATE(prw_ancien, prlw_ancien, prsw_ancien) -#ifdef ISO - REAL, ALLOCATABLE, SAVE :: xt_ancien(:,:,:),xtl_ancien(:,:,:),xts_ancien(:,:,:) -!$OMP THREADPRIVATE(xt_ancien,xtl_ancien,xts_ancien) -#endif - REAL, ALLOCATABLE, SAVE :: u_ancien(:,:), v_ancien(:,:) -!$OMP THREADPRIVATE(u_ancien, v_ancien) -!!! RomP >>> - REAL, ALLOCATABLE, SAVE :: tr_ancien(:,:,:) -!$OMP THREADPRIVATE(tr_ancien) -!!! RomP <<< - LOGICAL, SAVE :: ancien_ok -!$OMP THREADPRIVATE(ancien_ok) - REAL, ALLOCATABLE, SAVE :: clwcon(:,:),rnebcon(:,:) -!$OMP THREADPRIVATE(clwcon,rnebcon) - REAL, ALLOCATABLE, SAVE :: rneb_ancien(:,:) -!$OMP THREADPRIVATE(rneb_ancien) - REAL, ALLOCATABLE, SAVE :: qtc_cv(:,:),sigt_cv(:,:) -!$OMP THREADPRIVATE(qtc_cv,sigt_cv) - REAL, ALLOCATABLE, SAVE :: ratqs(:,:) -!$OMP THREADPRIVATE(ratqs) - REAL, ALLOCATABLE, SAVE :: pbl_tke(:,:,:) ! turb kinetic energy - REAL, ALLOCATABLE, SAVE :: coefh(:,:,:) ! Kz enthalpie - REAL, ALLOCATABLE, SAVE :: coefm(:,:,:) ! Kz momentum -!$OMP THREADPRIVATE(pbl_tke, coefh,coefm) - REAL, ALLOCATABLE, SAVE :: zmax0(:), f0(:) ! -!$OMP THREADPRIVATE(zmax0,f0) - REAL, ALLOCATABLE, SAVE :: sig1(:,:), w01(:,:) -!$OMP THREADPRIVATE(sig1,w01) - REAL, ALLOCATABLE, SAVE :: entr_therm(:,:), fm_therm(:,:) -!$OMP THREADPRIVATE(entr_therm,fm_therm) - REAL, ALLOCATABLE, SAVE :: detr_therm(:,:) -!$OMP THREADPRIVATE(detr_therm) -!IM 150408 -! pour phsystoke avec thermiques - REAL,ALLOCATABLE,SAVE :: clwcon0th(:,:),rnebcon0th(:,:) -!$OMP THREADPRIVATE(clwcon0th,rnebcon0th) -! radiation outputs - REAL,ALLOCATABLE,SAVE :: swdnc0(:,:), swdn0(:,:), swdn(:,:) -!$OMP THREADPRIVATE(swdnc0,swdn0,swdn) - REAL,ALLOCATABLE,SAVE :: swupc0(:,:), swup0(:,:), swup(:,:) -!$OMP THREADPRIVATE(swupc0, swup0,swup) - REAL,ALLOCATABLE,SAVE :: SWdn200clr(:), SWdn200(:) -!$OMP THREADPRIVATE(SWdn200clr,SWdn200) - REAL,ALLOCATABLE,SAVE :: SWup200clr(:), SWup200(:) -!$OMP THREADPRIVATE(SWup200clr,SWup200) - REAL,ALLOCATABLE,SAVE :: lwdnc0(:,:), lwdn0(:,:), lwdn(:,:) -!$OMP THREADPRIVATE(lwdnc0,lwdn0,lwdn) - REAL,ALLOCATABLE,SAVE :: lwupc0(:,:), lwup0(:,:), lwup(:,:) -!$OMP THREADPRIVATE(lwupc0,lwup0,lwup) - REAL,ALLOCATABLE,SAVE :: LWdn200clr(:), LWdn200(:) -!$OMP THREADPRIVATE(LWdn200clr,LWdn200) - REAL,ALLOCATABLE,SAVE :: LWup200clr(:), LWup200(:) -!$OMP THREADPRIVATE(LWup200clr,LWup200) - REAL,ALLOCATABLE,SAVE :: LWdnTOA(:), LWdnTOAclr(:) -!$OMP THREADPRIVATE(LWdnTOA,LWdnTOAclr) -! pressure level - REAL,ALLOCATABLE,SAVE :: tsumSTD(:,:,:) -!$OMP THREADPRIVATE(tsumSTD) - REAL,ALLOCATABLE,SAVE :: usumSTD(:,:,:), vsumSTD(:,:,:) -!$OMP THREADPRIVATE(usumSTD,vsumSTD) - REAL,ALLOCATABLE,SAVE :: wsumSTD(:,:,:), phisumSTD(:,:,:) -!$OMP THREADPRIVATE(wsumSTD,phisumSTD) - REAL,ALLOCATABLE,SAVE :: qsumSTD(:,:,:), rhsumSTD(:,:,:) -!$OMP THREADPRIVATE(qsumSTD,rhsumSTD) - REAL,ALLOCATABLE,SAVE :: tnondef(:,:,:) -!$OMP THREADPRIVATE(tnondef) - REAL,ALLOCATABLE,SAVE :: uvsumSTD(:,:,:) -!$OMP THREADPRIVATE(uvsumSTD) - REAL,ALLOCATABLE,SAVE :: vqsumSTD(:,:,:) -!$OMP THREADPRIVATE(vqsumSTD) - REAL,ALLOCATABLE,SAVE :: vTsumSTD(:,:,:) -!$OMP THREADPRIVATE(vTsumSTD) - REAL,ALLOCATABLE,SAVE :: wqsumSTD(:,:,:) -!$OMP THREADPRIVATE(wqsumSTD) - REAL,ALLOCATABLE,SAVE :: vphisumSTD(:,:,:) -!$OMP THREADPRIVATE(vphisumSTD) - REAL,ALLOCATABLE,SAVE :: wTsumSTD(:,:,:) -!$OMP THREADPRIVATE(wTsumSTD) - REAL,ALLOCATABLE,SAVE :: u2sumSTD(:,:,:) -!$OMP THREADPRIVATE(u2sumSTD) - REAL,ALLOCATABLE,SAVE :: v2sumSTD(:,:,:) -!$OMP THREADPRIVATE(v2sumSTD) - REAL,ALLOCATABLE,SAVE :: T2sumSTD(:,:,:) -!$OMP THREADPRIVATE(T2sumSTD) - REAL,ALLOCATABLE,SAVE :: O3sumSTD(:,:,:), O3daysumSTD(:,:,:) -!$OMP THREADPRIVATE(O3sumSTD,O3daysumSTD) -!IM begin - REAL,ALLOCATABLE,SAVE :: wlevSTD(:,:), ulevSTD(:,:), vlevSTD(:,:) -!$OMP THREADPRIVATE(wlevSTD,ulevSTD,vlevSTD) - REAL,ALLOCATABLE,SAVE :: tlevSTD(:,:), qlevSTD(:,:), rhlevSTD(:,:) -!$OMP THREADPRIVATE(tlevSTD,qlevSTD,rhlevSTD) - REAL,ALLOCATABLE,SAVE :: philevSTD(:,:) -!$OMP THREADPRIVATE(philevSTD) - REAL,ALLOCATABLE,SAVE :: uvSTD(:,:) -!$OMP THREADPRIVATE(uvSTD) - REAL,ALLOCATABLE,SAVE :: vqSTD(:,:) -!$OMP THREADPRIVATE(vqSTD) - REAL,ALLOCATABLE,SAVE :: vTSTD(:,:) -!$OMP THREADPRIVATE(vTSTD) - REAL,ALLOCATABLE,SAVE :: wqSTD(:,:) -!$OMP THREADPRIVATE(wqSTD) - REAL,ALLOCATABLE,SAVE :: vphiSTD(:,:) -!$OMP THREADPRIVATE(vphiSTD) - REAL,ALLOCATABLE,SAVE :: wTSTD(:,:) -!$OMP THREADPRIVATE(wTSTD) - REAL,ALLOCATABLE,SAVE :: u2STD(:,:) -!$OMP THREADPRIVATE(u2STD) - REAL,ALLOCATABLE,SAVE :: v2STD(:,:) -!$OMP THREADPRIVATE(v2STD) - REAL,ALLOCATABLE,SAVE :: T2STD(:,:) -!$OMP THREADPRIVATE(T2STD) - REAL,ALLOCATABLE,SAVE :: O3STD(:,:), O3daySTD(:,:) -!$OMP THREADPRIVATE(O3STD,O3daySTD) -!IM end - INTEGER,ALLOCATABLE,SAVE :: seed_old(:,:) -!$OMP THREADPRIVATE(seed_old) - REAL,ALLOCATABLE,SAVE :: zuthe(:),zvthe(:) -!$OMP THREADPRIVATE(zuthe,zvthe) - REAL,ALLOCATABLE,SAVE :: alb_neig(:) -!$OMP THREADPRIVATE(alb_neig) -!cloud base mass flux - REAL,ALLOCATABLE,SAVE :: ema_cbmf(:) -!$OMP THREADPRIVATE(ema_cbmf) -!cloud base pressure & cloud top pressure - REAL,ALLOCATABLE,SAVE :: ema_pcb(:), ema_pct(:) -!$OMP THREADPRIVATE(ema_pcb,ema_pct) - REAL,ALLOCATABLE,SAVE :: Mipsh(:,:) ! mass flux shed from adiab. ascents -!$OMP THREADPRIVATE(Mipsh) - REAL,ALLOCATABLE,SAVE :: Ma(:,:) ! undilute upward mass flux -!$OMP THREADPRIVATE(Ma) - REAL,ALLOCATABLE,SAVE :: qcondc(:,:) ! in-cld water content from convect -!$OMP THREADPRIVATE(qcondc) - REAL,ALLOCATABLE,SAVE :: wd(:) ! sb -!$OMP THREADPRIVATE(wd) - REAL,ALLOCATABLE,SAVE :: sigd(:) -!$OMP THREADPRIVATE(sigd) -! - REAL,ALLOCATABLE,SAVE :: cin(:) -!$OMP THREADPRIVATE(cin) -! ftd : convective heating due to unsaturated downdraughts - REAL,ALLOCATABLE,SAVE :: ftd(:,:) -!$OMP THREADPRIVATE(ftd) -! fqd : convective moistening due to unsaturated downdraughts - REAL,ALLOCATABLE,SAVE :: fqd(:,:) -!$OMP THREADPRIVATE(fqd) -#ifdef ISO - REAL, ALLOCATABLE, SAVE :: fxtd(:,:,:) -!$OMP THREADPRIVATE(fxtd) -#endif -!34EK -! -- Variables de controle de ALE et ALP -!ALE : Energie disponible pour soulevement : utilisee par la -! convection d'Emanuel pour le declenchement et la regulation - REAL,ALLOCATABLE,SAVE :: ALE(:) -!$OMP THREADPRIVATE(ALE) -!ALP : Puissance disponible pour soulevement - REAL,ALLOCATABLE,SAVE :: ALP(:) -!$OMP THREADPRIVATE(ALP) -! -! nouvelles variables pour le couplage convection-couche limite - REAL,ALLOCATABLE,SAVE :: Ale_bl(:) -!$OMP THREADPRIVATE(Ale_bl) - REAL,ALLOCATABLE,SAVE :: Alp_bl(:) -!$OMP THREADPRIVATE(Alp_bl) - INTEGER,ALLOCATABLE,SAVE :: lalim_conv(:) -!$OMP THREADPRIVATE(lalim_conv) - REAL,ALLOCATABLE,SAVE :: wght_th(:,:) -!$OMP THREADPRIVATE(wght_th) - REAL,ALLOCATABLE,SAVE :: ale_wake(:) -!$OMP THREADPRIVATE(ale_wake) - REAL,ALLOCATABLE,SAVE :: ale_bl_stat(:) -!$OMP THREADPRIVATE(ale_bl_stat) -! -! variables de la wake -! wake_deltat : ecart de temperature avec la zone non perturbee -! wake_deltaq : ecart d'humidite avec la zone non perturbee -! wake_s : fraction surfacique occupee par la poche froide -! awake_dens : number of active wakes per unit area -! wake_dens : number of wakes per unit area -! cv_gen : birth rate of cumulonimbus per unit area. -! wake_occ : occurence of wakes (= 1 if wakes occur, =0 otherwise) -! wake_Cstar : vitesse d'etalement de la poche -! wake_pe : wake potential energy - WAPE -! wake_fip : Gust Front Impinging power - ALP - REAL,ALLOCATABLE,SAVE :: wake_deltat(:,:) -!$OMP THREADPRIVATE(wake_deltat) - REAL,ALLOCATABLE,SAVE :: wake_deltaq(:,:) -!$OMP THREADPRIVATE(wake_deltaq) -#ifdef ISO - REAL, ALLOCATABLE, SAVE :: wake_deltaxt(:,:,:) -!$OMP THREADPRIVATE(wake_deltaxt) -#endif - REAL,ALLOCATABLE,SAVE :: wake_s(:) -!$OMP THREADPRIVATE(wake_s) - REAL,ALLOCATABLE,SAVE :: awake_dens(:), wake_dens(:) -!$OMP THREADPRIVATE(awake_dens, wake_dens) - REAL,ALLOCATABLE,SAVE :: cv_gen(:) -!$OMP THREADPRIVATE(cv_gen) - REAL,ALLOCATABLE,SAVE :: wake_Cstar(:) -!$OMP THREADPRIVATE(wake_Cstar) - REAL,ALLOCATABLE,SAVE :: wake_pe(:) -!$OMP THREADPRIVATE(wake_pe) - REAL,ALLOCATABLE,SAVE :: wake_fip(:) -!$OMP THREADPRIVATE(wake_fip) -! -!jyg< -! variables related to the spitting of the PBL between wake and -! off-wake regions. -! wake_delta_pbl_TKE : difference TKE_w - TKE_x - REAL,ALLOCATABLE,SAVE :: wake_delta_pbl_TKE(:,:,:) -!$OMP THREADPRIVATE(wake_delta_pbl_TKE) -!nrlmd< - REAL, ALLOCATABLE, SAVE :: delta_tsurf(:,:) ! Surface temperature difference inside-outside cold pool -!$OMP THREADPRIVATE(delta_tsurf) -!>nrlmd -!>jyg -! -! pfrac_impa : Produits des coefs lessivage impaction -! pfrac_nucl : Produits des coefs lessivage nucleation -! pfrac_1nucl: Produits des coefs lessi nucl (alpha = 1) - REAL,ALLOCATABLE,SAVE :: pfrac_impa(:,:), pfrac_nucl(:,:) -!$OMP THREADPRIVATE(pfrac_impa,pfrac_nucl) - REAL,ALLOCATABLE,SAVE :: pfrac_1nucl(:,:) -!$OMP THREADPRIVATE(pfrac_1nucl) -! - REAL,ALLOCATABLE,SAVE :: total_rain(:), nday_rain(:) -!$OMP THREADPRIVATE(total_rain,nday_rain) - REAL,ALLOCATABLE,SAVE :: paire_ter(:) -!$OMP THREADPRIVATE(paire_ter) -! albsol1: albedo du sol total pour SW visible -! albsol2: albedo du sol total pour SW proche IR - REAL,ALLOCATABLE,SAVE :: albsol1(:), albsol2(:) -!$OMP THREADPRIVATE(albsol1,albsol2) - -!albedo SB >>> - REAL,ALLOCATABLE,SAVE :: albsol_dif(:,:),albsol_dir(:,:) -!$OMP THREADPRIVATE(albsol_dif,albsol_dir) -!albedo SB <<< - - - REAL, ALLOCATABLE, SAVE:: wo(:, :, :) - ! column-density of ozone in a layer, in kilo-Dobsons - ! Third dimension has size 1 or 2. - ! "wo(:, :, 1)" is for the average day-night field, - ! "wo(:, :, 2)" is for daylight time. - !$OMP THREADPRIVATE(wo) - -! heat : chauffage solaire -! heat0: chauffage solaire ciel clair -! cool : refroidissement infrarouge -! cool0 : refroidissement infrarouge ciel clair -! sollwdown : downward LW flux at surface -! sollwdownclr : downward CS LW flux at surface -! toplwdown : downward CS LW flux at TOA -! toplwdownclr : downward CS LW flux at TOA -! heat_volc : chauffage solaire du au volcanisme -! cool_volc : refroidissement infrarouge du au volcanisme - REAL,ALLOCATABLE,SAVE :: clwcon0(:,:),rnebcon0(:,:) -!$OMP THREADPRIVATE(clwcon0,rnebcon0) - REAL,ALLOCATABLE,SAVE :: heat(:,:) -!$OMP THREADPRIVATE(heat) - REAL,ALLOCATABLE,SAVE :: heat0(:,:) -!$OMP THREADPRIVATE(heat0) - REAL,ALLOCATABLE,SAVE :: cool(:,:) -!$OMP THREADPRIVATE(cool) - REAL,ALLOCATABLE,SAVE :: cool0(:,:) -!$OMP THREADPRIVATE(cool0) - REAL,ALLOCATABLE,SAVE :: heat_volc(:,:) -!$OMP THREADPRIVATE(heat_volc) - REAL,ALLOCATABLE,SAVE :: cool_volc(:,:) -!$OMP THREADPRIVATE(cool_volc) - REAL,ALLOCATABLE,SAVE :: topsw(:), toplw(:) -!$OMP THREADPRIVATE(topsw,toplw) - REAL,ALLOCATABLE,SAVE :: sollwdown(:) -!$OMP THREADPRIVATE(sollwdown) - REAL,ALLOCATABLE,SAVE :: gustiness(:) -!$OMP THREADPRIVATE(gustiness) - REAL,ALLOCATABLE,SAVE :: sollwdownclr(:) -!$OMP THREADPRIVATE(sollwdownclr) - REAL,ALLOCATABLE,SAVE :: toplwdown(:) -!$OMP THREADPRIVATE(toplwdown) - REAL,ALLOCATABLE,SAVE :: toplwdownclr(:) -!$OMP THREADPRIVATE(toplwdownclr) - REAL,ALLOCATABLE,SAVE :: topsw0(:),toplw0(:),solsw0(:),sollw0(:) -!$OMP THREADPRIVATE(topsw0,toplw0,solsw0,sollw0) - REAL,ALLOCATABLE,SAVE :: albpla(:) -!$OMP THREADPRIVATE(albpla) - -!IM ajout variables CFMIP2/CMIP5 - REAL,ALLOCATABLE,SAVE :: heatp(:,:), coolp(:,:) -!$OMP THREADPRIVATE(heatp, coolp) - REAL,ALLOCATABLE,SAVE :: heat0p(:,:), cool0p(:,:) -!$OMP THREADPRIVATE(heat0p, cool0p) - REAL,ALLOCATABLE,SAVE :: radsolp(:), topswp(:), toplwp(:) -!$OMP THREADPRIVATE(radsolp, topswp, toplwp) - REAL,ALLOCATABLE,SAVE :: albplap(:) -!$OMP THREADPRIVATE(albplap) - REAL,ALLOCATABLE,SAVE :: solswp(:), solswfdiffp(:), sollwp(:) -!$OMP THREADPRIVATE(solswp, solswfdiffp, sollwp) - REAL,ALLOCATABLE,SAVE :: sollwdownp(:) -!$OMP THREADPRIVATE(sollwdownp) - REAL,ALLOCATABLE,SAVE :: topsw0p(:),toplw0p(:) - REAL,ALLOCATABLE,SAVE :: solsw0p(:),sollw0p(:) -!$OMP THREADPRIVATE(topsw0p,toplw0p,solsw0p,sollw0p) - REAL,ALLOCATABLE,SAVE :: lwdnc0p(:,:), lwdn0p(:,:), lwdnp(:,:) - REAL,ALLOCATABLE,SAVE :: lwupc0p(:,:), lwup0p(:,:), lwupp(:,:) -!$OMP THREADPRIVATE(lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp) - REAL,ALLOCATABLE,SAVE :: swdnc0p(:,:), swdn0p(:,:), swdnp(:,:) - REAL,ALLOCATABLE,SAVE :: swupc0p(:,:), swup0p(:,:), swupp(:,:) -!$OMP THREADPRIVATE(swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp) - -! pbase : cloud base pressure -! bbase : cloud base buoyancy - REAL,ALLOCATABLE,SAVE :: cape(:) -!$OMP THREADPRIVATE(cape) - REAL,ALLOCATABLE,SAVE :: pbase(:) -!$OMP THREADPRIVATE(pbase) - REAL,ALLOCATABLE,SAVE :: bbase(:) -!$OMP THREADPRIVATE(bbase) -! - REAL,SAVE,ALLOCATABLE :: zqasc(:,:) -!$OMP THREADPRIVATE( zqasc) - INTEGER,ALLOCATABLE,SAVE :: ibas_con(:), itop_con(:) -!$OMP THREADPRIVATE(ibas_con,itop_con) - REAL,SAVE,ALLOCATABLE :: rain_con(:) -!$OMP THREADPRIVATE(rain_con) - REAL,SAVE,ALLOCATABLE :: snow_con(:) -!$OMP THREADPRIVATE(snow_con) -! -#ifdef ISO - REAL,SAVE,ALLOCATABLE :: xtrain_con(:,:) -!$OMP THREADPRIVATE(xtrain_con) - REAL,SAVE,ALLOCATABLE :: xtsnow_con(:,:) -!$OMP THREADPRIVATE(xtsnow_con) -#endif - REAL,SAVE,ALLOCATABLE :: rlonPOS(:) -!$OMP THREADPRIVATE(rlonPOS) - REAL,SAVE,ALLOCATABLE :: newsst(:) -!$OMP THREADPRIVATE(newsst) - REAL,SAVE,ALLOCATABLE :: ustar(:,:),u10m(:,:), v10m(:,:),wstar(:,:) -!$OMP THREADPRIVATE(ustar,u10m,v10m,wstar) -! -! ok_ade=T -ADE=topswad-topsw -! ok_aie=T -> -! ok_ade=T -AIE=topswai-topswad -! ok_ade=F -AIE=topswai-topsw -! -!topswad, solswad : Aerosol direct effect - REAL,SAVE,ALLOCATABLE :: topswad(:), solswad(:) -!$OMP THREADPRIVATE(topswad,solswad) -!topswai, solswai : Aerosol indirect effect - REAL,SAVE,ALLOCATABLE :: topswai(:), solswai(:) -!$OMP THREADPRIVATE(topswai,solswai) - - REAL,SAVE,ALLOCATABLE :: tau_aero(:,:,:,:), piz_aero(:,:,:,:), cg_aero(:,:,:,:) -!$OMP THREADPRIVATE(tau_aero, piz_aero, cg_aero) - REAL,SAVE,ALLOCATABLE :: tau_aero_sw_rrtm(:,:,:,:), piz_aero_sw_rrtm(:,:,:,:), cg_aero_sw_rrtm(:,:,:,:) -!$OMP THREADPRIVATE(tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm) - REAL,SAVE,ALLOCATABLE :: tau_aero_lw_rrtm(:,:,:,:), piz_aero_lw_rrtm(:,:,:,:), cg_aero_lw_rrtm(:,:,:,:) -!$OMP THREADPRIVATE(tau_aero_lw_rrtm, piz_aero_lw_rrtm, cg_aero_lw_rrtm) - REAL,SAVE,ALLOCATABLE :: ccm(:,:,:) -!$OMP THREADPRIVATE(ccm) - - REAL,SAVE,ALLOCATABLE :: ale_bl_trig(:) -!$OMP THREADPRIVATE(ale_bl_trig) - - REAL,SAVE,ALLOCATABLE :: ratqs_inter(:,:) -!$OMP THREADPRIVATE(ratqs_inter) - -#ifdef ISO -#ifdef ISOTRAC - INTEGER,SAVE,ALLOCATABLE :: bassin_map(:) -!$OMP THREADPRIVATE(bassin_map) - INTEGER,SAVE,ALLOCATABLE :: boite_map(:,:) -!$OMP THREADPRIVATE(boite_map) -#endif -#endif - REAL, ALLOCATABLE, SAVE:: du_gwd_rando(:, :), du_gwd_front(:, :) - !$OMP THREADPRIVATE(du_gwd_rando, du_gwd_front) - ! tendencies on wind due to gravity waves - - LOGICAL,SAVE :: is_initialized=.FALSE. -!$OMP THREADPRIVATE(is_initialized) - - ! Ocean-atmosphere interface: - - REAL, ALLOCATABLE, SAVE:: ds_ns(:) ! (klon) - ! "delta salinity near surface". Salinity variation in the - ! near-surface turbulent layer. That is subskin salinity minus - ! foundation salinity. In ppt. - - REAL, ALLOCATABLE, SAVE:: dt_ns(:) ! (klon) - ! "delta temperature near surface". Temperature variation in the - ! near-surface turbulent layer. That is subskin temperature - ! minus foundation temperature. (Can be negative.) In K. - - REAL, ALLOCATABLE, SAVE:: delta_sst(:) ! (klon) - ! Ocean-air interface temperature minus bulk SST, in - ! K. Allocated and defined only if activate_ocean_skin >= 1. - - REAL, ALLOCATABLE, SAVE:: delta_sal(:) ! (klon) - ! Ocean-air interface salinity minus bulk salinity, in ppt - - !$OMP THREADPRIVATE(delta_sal, ds_ns, dt_ns, delta_sst) - - CONTAINS - -!====================================================================== -SUBROUTINE phys_state_var_init(read_climoz) -USE dimphy -USE aero_mod -USE infotrac_phy, ONLY : nbtr -#ifdef ISO -USE infotrac_phy, ONLY : ntraciso,niso -#endif -USE indice_sol_mod -use config_ocean_skin_m, only: activate_ocean_skin -IMPLICIT NONE - -integer, intent(in):: read_climoz -! read ozone climatology -! Allowed values are 0, 1 and 2 -! 0: do not read an ozone climatology -! 1: read a single ozone climatology that will be used day and night -! 2: read two ozone climatologies, the average day and night -! climatology and the daylight climatology - -include "clesphys.h" - - print*, 'is_initialized', is_initialized - IF (is_initialized) RETURN - is_initialized=.TRUE. - ALLOCATE(pctsrf(klon,nbsrf)) - ALLOCATE(ftsol(klon,nbsrf)) - ALLOCATE(beta_aridity(klon,nbsrf)) - ALLOCATE(qsol(klon),fevap(klon,nbsrf)) - ALLOCATE(z0m(klon,nbsrf+1),z0h(klon,nbsrf+1),agesno(klon,nbsrf)) -!FC - ALLOCATE(treedrg(klon,klev,nbsrf)) - ALLOCATE(falb1(klon,nbsrf)) - ALLOCATE(falb2(klon,nbsrf)) -!albedo SB >>> - print*, 'allocate falb' - ALLOCATE(falb_dir(klon,nsw,nbsrf),falb_dif(klon,nsw,nbsrf)) -!! print*, 'allocate falb good', falb_dir(1,1,1) - ALLOCATE(chl_con(klon)) -!albedo SB <<< - ALLOCATE(rain_fall(klon)) - ALLOCATE(snow_fall(klon)) - ALLOCATE(solsw(klon), solswfdiff(klon), sollw(klon)) - sollw=0.0 - ALLOCATE(radsol(klon)) - ALLOCATE(swradcorr(klon)) - ALLOCATE(zmea(klon), zstd(klon), zsig(klon), zgam(klon)) - ALLOCATE(zthe(klon), zpic(klon), zval(klon)) - - ALLOCATE(rugoro(klon)) - ALLOCATE(t_ancien(klon,klev), q_ancien(klon,klev)) - ALLOCATE(ql_ancien(klon,klev), qs_ancien(klon,klev)) - ALLOCATE(prw_ancien(klon), prlw_ancien(klon), prsw_ancien(klon)) - ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev)) -!!! Rom P >>> - ALLOCATE(tr_ancien(klon,klev,nbtr)) -!!! Rom P <<< - ALLOCATE(clwcon(klon,klev),rnebcon(klon,klev)) - ALLOCATE(rneb_ancien(klon,klev)) - ALLOCATE(qtc_cv(klon,klev),sigt_cv(klon,klev)) - ALLOCATE(ratqs(klon,klev)) - ALLOCATE(pbl_tke(klon,klev+1,nbsrf+1)) -!nrlmd< - ALLOCATE(delta_tsurf(klon,nbsrf)) -!>nrlmd - ALLOCATE(coefh(klon,klev+1,nbsrf+1)) - ALLOCATE(coefm(klon,klev+1,nbsrf+1)) - ALLOCATE(zmax0(klon), f0(klon)) - ALLOCATE(sig1(klon,klev), w01(klon,klev)) - ALLOCATE(entr_therm(klon,klev), fm_therm(klon,klev+1)) - ALLOCATE(detr_therm(klon,klev)) -! pour phsystoke avec thermiques - ALLOCATE(clwcon0th(klon,klev),rnebcon0th(klon,klev)) -! radiation outputs - ALLOCATE(swdnc0(klon,klevp1), swdn0(klon,klevp1), swdn(klon,klevp1)) - ALLOCATE(swupc0(klon,klevp1), swup0(klon,klevp1), swup(klon,klevp1)) - ALLOCATE(lwdnc0(klon,klevp1), lwdn0(klon,klevp1), lwdn(klon,klevp1)) - ALLOCATE(lwupc0(klon,klevp1), lwup0(klon,klevp1), lwup(klon,klevp1)) - ALLOCATE(SWdn200clr(klon), SWdn200(klon)) - ALLOCATE(SWup200clr(klon), SWup200(klon)) - ALLOCATE(LWdn200clr(klon), LWdn200(klon)) - ALLOCATE(LWup200clr(klon), LWup200(klon)) - ALLOCATE(LWdnTOA(klon), LWdnTOAclr(klon)) -! pressure level - ALLOCATE(tsumSTD(klon,nlevSTD,nout)) - ALLOCATE(usumSTD(klon,nlevSTD,nout), vsumSTD(klon,nlevSTD,nout)) - ALLOCATE(wsumSTD(klon,nlevSTD,nout), phisumSTD(klon,nlevSTD,nout)) - ALLOCATE(qsumSTD(klon,nlevSTD,nout), rhsumSTD(klon,nlevSTD,nout)) - ALLOCATE(tnondef(klon,nlevSTD,nout)) - ALLOCATE(uvsumSTD(klon,nlevSTD,nout)) - ALLOCATE(vqsumSTD(klon,nlevSTD,nout)) - ALLOCATE(vTsumSTD(klon,nlevSTD,nout)) - ALLOCATE(wqsumSTD(klon,nlevSTD,nout)) - ALLOCATE(vphisumSTD(klon,nlevSTD,nout)) - ALLOCATE(wTsumSTD(klon,nlevSTD,nout)) - ALLOCATE(u2sumSTD(klon,nlevSTD,nout)) - ALLOCATE(v2sumSTD(klon,nlevSTD,nout)) - ALLOCATE(T2sumSTD(klon,nlevSTD,nout)) - ALLOCATE(O3sumSTD(klon,nlevSTD,nout)) - ALLOCATE(O3daysumSTD(klon,nlevSTD,nout)) -!IM beg - ALLOCATE(wlevSTD(klon,nlevSTD), ulevSTD(klon,nlevSTD), vlevSTD(klon,nlevSTD)) - ALLOCATE(tlevSTD(klon,nlevSTD), qlevSTD(klon,nlevSTD), rhlevSTD(klon,nlevSTD)) - ALLOCATE(philevSTD(klon,nlevSTD)) - ALLOCATE(uvSTD(klon,nlevSTD),vqSTD(klon,nlevSTD)) - ALLOCATE(vTSTD(klon,nlevSTD),wqSTD(klon,nlevSTD)) - ALLOCATE(vphiSTD(klon,nlevSTD),wTSTD(klon,nlevSTD)) - ALLOCATE(u2STD(klon,nlevSTD),v2STD(klon,nlevSTD)) - ALLOCATE(T2STD(klon,nlevSTD)) - ALLOCATE(O3STD(klon,nlevSTD)) - ALLOCATE(O3daySTD(klon,nlevSTD)) -!IM end - ALLOCATE(seed_old(klon,napisccp)) - ALLOCATE(zuthe(klon),zvthe(klon)) - ALLOCATE(alb_neig(klon)) -!cloud base mass flux - ALLOCATE(ema_cbmf(klon)) -!cloud base pressure & cloud top pressure - ALLOCATE(ema_pcb(klon), ema_pct(klon)) -! - ALLOCATE(Mipsh(klon,klev)) - ALLOCATE(Ma(klon,klev)) - ALLOCATE(qcondc(klon,klev)) - ALLOCATE(wd(klon)) - ALLOCATE(sigd(klon)) - ALLOCATE(cin(klon), ALE(klon), ALP(klon)) - ALLOCATE(ftd(klon,klev), fqd(klon,klev)) - ALLOCATE(Ale_bl(klon)) - ALLOCATE(ale_wake(klon)) - ALLOCATE(ale_bl_stat(klon)) - ale_bl_stat(:)=0 - ALLOCATE(Alp_bl(klon)) - ALLOCATE(lalim_conv(klon)) - ALLOCATE(wght_th(klon,klev)) - ALLOCATE(wake_deltat(klon,klev), wake_deltaq(klon,klev)) - ALLOCATE(wake_s(klon), awake_dens(klon), wake_dens(klon)) -!! awake_dens = 0. ! initialized in phyetat0 - ALLOCATE(cv_gen(klon)) - ALLOCATE(wake_Cstar(klon)) - ALLOCATE(wake_pe(klon), wake_fip(klon)) -!jyg< - ALLOCATE(wake_delta_pbl_TKE(klon,klev+1,nbsrf+1)) -!>jyg - ALLOCATE(pfrac_impa(klon,klev), pfrac_nucl(klon,klev)) - ALLOCATE(pfrac_1nucl(klon,klev)) - ALLOCATE(total_rain(klon), nday_rain(klon)) - ALLOCATE(paire_ter(klon)) - ALLOCATE(albsol1(klon), albsol2(klon)) -!albedo SB >>> - ALLOCATE(albsol_dir(klon,nsw),albsol_dif(klon,nsw)) -!albedo SB <<< - - if (read_climoz <= 1) then - ALLOCATE(wo(klon,klev, 1)) - else - ! read_climoz == 2 - ALLOCATE(wo(klon,klev, 2)) - end if - - ALLOCATE(clwcon0(klon,klev),rnebcon0(klon,klev)) - ALLOCATE(heat(klon,klev), heat0(klon,klev)) - ALLOCATE(cool(klon,klev), cool0(klon,klev)) - ALLOCATE(heat_volc(klon,klev), cool_volc(klon,klev)) - ALLOCATE(topsw(klon), toplw(klon)) - ALLOCATE(sollwdown(klon), sollwdownclr(klon)) - sollwdown = 0. - ALLOCATE(toplwdown(klon), toplwdownclr(klon)) - ALLOCATE(topsw0(klon),toplw0(klon),solsw0(klon),sollw0(klon)) - sollw0 = 0. - ALLOCATE(albpla(klon)) -!IM ajout variables CFMIP2/CMIP5 - ALLOCATE(heatp(klon,klev), coolp(klon,klev)) - ALLOCATE(heat0p(klon,klev), cool0p(klon,klev)) - ALLOCATE(radsolp(klon), topswp(klon), toplwp(klon)) - ALLOCATE(albplap(klon)) - ALLOCATE(solswp(klon), solswfdiffp(klon), sollwp(klon)) - ALLOCATE(gustiness(klon)) - ALLOCATE(sollwdownp(klon)) - ALLOCATE(topsw0p(klon),toplw0p(klon)) - ALLOCATE(solsw0p(klon),sollw0p(klon)) - ALLOCATE(lwdnc0p(klon,klevp1), lwdn0p(klon,klevp1), lwdnp(klon,klevp1)) - ALLOCATE(lwupc0p(klon,klevp1), lwup0p(klon,klevp1), lwupp(klon,klevp1)) - ALLOCATE(swdnc0p(klon,klevp1), swdn0p(klon,klevp1), swdnp(klon,klevp1)) - ALLOCATE(swupc0p(klon,klevp1), swup0p(klon,klevp1), swupp(klon,klevp1)) - - ALLOCATE(cape(klon)) - ALLOCATE(pbase(klon),bbase(klon)) - ALLOCATE(zqasc(klon,klev)) - ALLOCATE(ibas_con(klon), itop_con(klon)) - ALLOCATE(rain_con(klon), snow_con(klon)) - ALLOCATE(rlonPOS(klon)) - ALLOCATE(newsst(klon)) - ALLOCATE(ustar(klon,nbsrf),u10m(klon,nbsrf), v10m(klon,nbsrf),wstar(klon,nbsrf+1)) - ALLOCATE(topswad(klon), solswad(klon)) - ALLOCATE(topswai(klon), solswai(klon)) - ALLOCATE(tau_aero(klon,klev,naero_grp,nbands),piz_aero(klon,klev,naero_grp,nbands),cg_aero(klon,klev,naero_grp,nbands)) - ALLOCATE(tau_aero_sw_rrtm(klon,klev,2,nbands_sw_rrtm),piz_aero_sw_rrtm(klon,klev,2,nbands_sw_rrtm)) - ALLOCATE(cg_aero_sw_rrtm(klon,klev,2,nbands_sw_rrtm)) - ALLOCATE(tau_aero_lw_rrtm(klon,klev,2,nbands_lw_rrtm),piz_aero_lw_rrtm(klon,klev,2,nbands_lw_rrtm)) - ALLOCATE(cg_aero_lw_rrtm(klon,klev,2,nbands_lw_rrtm)) - ALLOCATE(ccm(klon,klev,nbands)) - -#ifdef ISO - ALLOCATE(xtsol(niso,klon),fxtevap(ntraciso,klon,nbsrf)) - ALLOCATE(fxtd(ntraciso,klon,klev)) - ALLOCATE(wake_deltaxt(ntraciso,klon,klev)) - ALLOCATE(xt_ancien(ntraciso,klon,klev)) - ALLOCATE(xtl_ancien(ntraciso,klon,klev)) - ALLOCATE(xts_ancien(ntraciso,klon,klev)) - ALLOCATE(xtrain_fall(ntraciso,klon)) - ALLOCATE(xtsnow_fall(ntraciso,klon)) - ALLOCATE(xtrain_con(ntraciso,klon)) - ALLOCATE(xtsnow_con(ntraciso,klon)) -#ifdef ISOTRAC - ALLOCATE(bassin_map(klon)) - ALLOCATE(boite_map(klon,klev)) -#endif -#endif - - ALLOCATE(ale_bl_trig(klon)) - ALLOCATE(ratqs_inter(klon,klev)) - IF (ok_gwd_rando) THEN - ALLOCATE(du_gwd_rando(klon, klev)) - du_gwd_rando(:,:)=0. - ENDIF - IF (.not. ok_hines .and. ok_gwd_rando) THEN - ALLOCATE(du_gwd_front(klon, klev)) - du_gwd_front(:,:) = 0 !ym missing init - ENDIF - if (activate_ocean_skin >= 1) ALLOCATE(delta_sal(klon), ds_ns(klon), & - dt_ns(klon), delta_sst(klon)) - - END SUBROUTINE phys_state_var_init - -!====================================================================== - SUBROUTINE phys_state_var_end - ! Useful only for lmdz1d. -!USE dimphy -USE indice_sol_mod -use config_ocean_skin_m, only: activate_ocean_skin -IMPLICIT NONE -include "clesphys.h" - - DEALLOCATE(pctsrf, ftsol, falb1, falb2) - DEALLOCATE(beta_aridity) - DEALLOCATE(qsol,fevap,z0m,z0h,agesno) -!FC - DEALLOCATE(treedrg) - DEALLOCATE(rain_fall, snow_fall, solsw, solswfdiff, sollw, radsol, swradcorr) - DEALLOCATE(zmea, zstd, zsig, zgam) - DEALLOCATE(zthe, zpic, zval) - DEALLOCATE(rugoro, t_ancien, q_ancien, clwcon, rnebcon) - DEALLOCATE(qs_ancien, ql_ancien, rneb_ancien) - DEALLOCATE(prw_ancien, prlw_ancien, prsw_ancien) - DEALLOCATE(qtc_cv,sigt_cv) - DEALLOCATE(u_ancien, v_ancien) - DEALLOCATE(tr_ancien) !RomP - DEALLOCATE(ratqs, pbl_tke,coefh,coefm) - DEALLOCATE(zmax0, f0) - DEALLOCATE(sig1, w01) - DEALLOCATE(entr_therm, fm_therm) - DEALLOCATE(detr_therm) - DEALLOCATE(clwcon0th, rnebcon0th) -! radiation outputs - DEALLOCATE(swdnc0, swdn0, swdn) - DEALLOCATE(swupc0, swup0, swup) - DEALLOCATE(lwdnc0, lwdn0, lwdn) - DEALLOCATE(lwupc0, lwup0, lwup) - DEALLOCATE(SWdn200clr, SWdn200) - DEALLOCATE(SWup200clr, SWup200) - DEALLOCATE(LWdn200clr, LWdn200) - DEALLOCATE(LWup200clr, LWup200) - DEALLOCATE(LWdnTOA, LWdnTOAclr) -! pressure level - DEALLOCATE(tsumSTD) - DEALLOCATE(usumSTD, vsumSTD) - DEALLOCATE(wsumSTD, phisumSTD) - DEALLOCATE(tnondef) - DEALLOCATE(qsumSTD, rhsumSTD) - DEALLOCATE(uvsumSTD) - DEALLOCATE(vqsumSTD) - DEALLOCATE(vTsumSTD) - DEALLOCATE(wqsumSTD) - DEALLOCATE(vphisumSTD) - DEALLOCATE(wTsumSTD) - DEALLOCATE(u2sumSTD) - DEALLOCATE(v2sumSTD) - DEALLOCATE(T2sumSTD) - DEALLOCATE(O3sumSTD) - DEALLOCATE(O3daysumSTD) -!IM beg - DEALLOCATE(wlevSTD,ulevSTD,vlevSTD,tlevSTD,qlevSTD,rhlevSTD,philevSTD) - DEALLOCATE(uvSTD,vqSTD,vTSTD,wqSTD,vphiSTD,wTSTD,u2STD,v2STD,T2STD,O3STD,O3daySTD) -!IM end - DEALLOCATE(seed_old) - DEALLOCATE(zuthe, zvthe) - DEALLOCATE(alb_neig) - DEALLOCATE(ema_cbmf) - DEALLOCATE(ema_pcb, ema_pct) - DEALLOCATE(Mipsh, Ma, qcondc) - DEALLOCATE(wd, sigd) - DEALLOCATE(cin, ALE, ALP) - DEALLOCATE(ftd, fqd) - DEALLOCATE(Ale_bl, Alp_bl) - DEALLOCATE(ale_wake) - DEALLOCATE(ale_bl_stat) - DEALLOCATE(lalim_conv, wght_th) - DEALLOCATE(wake_deltat, wake_deltaq) - DEALLOCATE(wake_s, awake_dens, wake_dens) - DEALLOCATE(cv_gen) - DEALLOCATE(wake_Cstar, wake_pe, wake_fip) -!jyg< - DEALLOCATE(wake_delta_pbl_TKE) -!nrlmd< - DEALLOCATE(delta_tsurf) -!>nrlmd -!>jyg - DEALLOCATE(pfrac_impa, pfrac_nucl) - DEALLOCATE(pfrac_1nucl) - DEALLOCATE(total_rain, nday_rain) - DEALLOCATE(paire_ter) - DEALLOCATE(albsol1, albsol2) -!albedo SB >>> - DEALLOCATE(albsol_dir,albsol_dif,falb_dir,falb_dif,chl_con) -!albedo SB <<< - DEALLOCATE(wo) - DEALLOCATE(clwcon0,rnebcon0) - DEALLOCATE(heat, heat0) - DEALLOCATE(cool, cool0) - DEALLOCATE(heat_volc, cool_volc) - DEALLOCATE(topsw, toplw) - DEALLOCATE(sollwdown, sollwdownclr) - DEALLOCATE(gustiness) - DEALLOCATE(toplwdown, toplwdownclr) - DEALLOCATE(topsw0,toplw0,solsw0,sollw0) - DEALLOCATE(albpla) -!IM ajout variables CFMIP2/CMIP5 - DEALLOCATE(heatp, coolp) - DEALLOCATE(heat0p, cool0p) - DEALLOCATE(radsolp, topswp, toplwp) - DEALLOCATE(albplap) - DEALLOCATE(solswp, solswfdiffp, sollwp) - DEALLOCATE(sollwdownp) - DEALLOCATE(topsw0p,toplw0p) - DEALLOCATE(solsw0p,sollw0p) - DEALLOCATE(lwdnc0p, lwdn0p, lwdnp) - DEALLOCATE(lwupc0p, lwup0p, lwupp) - DEALLOCATE(swdnc0p, swdn0p, swdnp) - DEALLOCATE(swupc0p, swup0p, swupp) - DEALLOCATE(cape) - DEALLOCATE(pbase,bbase) - DEALLOCATE(zqasc) - DEALLOCATE(ibas_con, itop_con) - DEALLOCATE(rain_con, snow_con) - DEALLOCATE(rlonPOS) - DEALLOCATE(newsst) - DEALLOCATE(ustar,u10m, v10m,wstar) - DEALLOCATE(topswad, solswad) - DEALLOCATE(topswai, solswai) - DEALLOCATE(tau_aero,piz_aero,cg_aero) - DEALLOCATE(tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm) - DEALLOCATE(tau_aero_lw_rrtm,piz_aero_lw_rrtm,cg_aero_lw_rrtm) - DEALLOCATE(ccm) - if (ok_gwd_rando) DEALLOCATE(du_gwd_rando) - if (.not. ok_hines .and. ok_gwd_rando) DEALLOCATE(du_gwd_front) - DEALLOCATE(ale_bl_trig) - DEALLOCATE(ratqs_inter) - - if (activate_ocean_skin >= 1) deALLOCATE(delta_sal, ds_ns, dt_ns, & - delta_sst) - -#ifdef ISO - DEALLOCATE(xtsol,fxtevap) - DEALLOCATE(xt_ancien,xtl_ancien,xts_ancien, fxtd, wake_deltaxt) - DEALLOCATE(xtrain_fall, xtsnow_fall, xtrain_con, xtsnow_con) -#ifdef ISOTRAC - DEALLOCATE(bassin_map,boite_map) -#endif -#endif - is_initialized=.FALSE. - -END SUBROUTINE phys_state_var_end - - END MODULE phys_state_var_mod Index: LMDZ6/trunk/libf/phylmdiso/thermcell_main.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/thermcell_main.F90 (revision 4090) +++ (revision ) @@ -1,857 +1,0 @@ - -! $Id: thermcell_main.F90 3451 2019-01-27 11:07:30Z fhourdin $ -! - 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 & - & ,zmax0, f0,zw2,fraca,ztv & - & ,zpspsk,ztla,zthl,ztva & - & ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax & -#ifdef ISO - & ,xtpo,xtpdoadj & -#endif - & ) - - - USE thermcell_ini_mod, ONLY: thermcell_ini,dqimpl,dvdq,prt_level,lunout,prt_level - USE thermcell_ini_mod, ONLY: iflag_thermals_closure,iflag_thermals_ed,tau_thermals,r_aspect_thermals - USE thermcell_ini_mod, ONLY: RD,RG - -#ifdef ISO - USE infotrac_phy, ONLY : ntraciso -#ifdef ISOVERIF - USE isotopes_mod, ONLY : iso_eau,iso_HDO - USE isotopes_verif_mod, ONLY: iso_verif_egalite, & - iso_verif_aberrant_encadre -#endif -#endif - - - 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: -! ------------- - - -! arguments: -! ---------- - integer, intent(in) :: itap,ngrid,nlay - real, intent(in) :: ptimestep - real, intent(in), dimension(ngrid,nlay) :: pt,pu,pv,po,pplay,pphi,zpspsk - real, intent(in), dimension(ngrid,nlay+1) :: pplev - real, intent(out), dimension(ngrid,nlay) :: pdtadj,pduadj,pdvadj,pdoadj,entr0,detr0 - real, intent(out), dimension(ngrid,nlay) :: ztla,zqla,zqta,zqsatth,zthl - real, intent(out), dimension(ngrid,nlay+1) :: fm0,zw2,fraca - real, intent(out), dimension(ngrid) :: zmax0,f0 - real, intent(out), dimension(ngrid,nlay) :: ztva,ztv - logical, intent(in) :: debut - - real, intent(out), dimension(ngrid) :: pcon - real, intent(out), dimension(ngrid,nlay) :: rhobarz,wth3 - real, intent(out), dimension(ngrid) :: wmax_sec - integer,intent(out), dimension(ngrid) :: lalim - real, intent(out), dimension(ngrid,nlay+1) :: fm - real, intent(out), dimension(ngrid,nlay) :: alim_star - real, intent(out), dimension(ngrid) :: zmax - -! local: -! ------ - - - integer,save :: igout=1 -!$OMP THREADPRIVATE(igout) - integer,save :: lunout1=6 -!$OMP THREADPRIVATE(lunout1) - integer,save :: lev_out=10 -!$OMP THREADPRIVATE(lev_out) - - real lambda, zf,zf2,var,vardiff,CHI - integer ig,k,l,ierr,ll - logical sorties - real, dimension(ngrid) :: linter,zmix, zmax_sec - integer,dimension(ngrid) :: lmax,lmin,lmix,lmix_bis,nivcon - real, dimension(ngrid,nlay) :: ztva_est - real, dimension(ngrid,nlay) :: deltaz,zlay,zh,zdthladj,zu,zv,zo,zl,zva,zua,zoa - real, dimension(ngrid,nlay) :: zta,zha,q2,wq,wthl,wthv,thetath2,wth2 - real, dimension(ngrid,nlay) :: ratqscth,ratqsdiff,rho,masse - real, dimension(ngrid,nlay+1) :: zw_est,zlev - real, dimension(ngrid) :: wmax,wmax_tmp - real, dimension(ngrid,nlay+1) :: f_star - real, dimension(ngrid,nlay) :: entr,detr,entr_star,detr_star,alim_star_clos - real, dimension(ngrid,nlay) :: zqsat,csc - real, dimension(ngrid) :: zcon,zcon2,alim_star_tot,f - - character (len=20) :: modname='thermcell_main' - character (len=80) :: abort_message - - -#ifdef ISO - REAL xtpo(ntraciso,ngrid,nlay),xtpdoadj(ntraciso,ngrid,nlay) - REAL xtzo(ntraciso,ngrid,nlay) - REAL xtpdoadj_tmp(ngrid,nlay) - REAL xtpo_tmp(ngrid,nlay) - REAL xtzo_tmp(ngrid,nlay) - integer ixt -#endif - -! - -!----------------------------------------------------------------------- -! initialisation: -! --------------- -! - print*,'NEW THERMCELL cool' - - - fm=0. ; entr=0. ; detr=0. - - if (prt_level.ge.1) print*,'thermcell_main V4' - - sorties=.true. - IF(ngrid.NE.ngrid) THEN - PRINT* - PRINT*,'STOP dans convadj' - PRINT*,'ngrid =',ngrid - PRINT*,'ngrid =',ngrid - ENDIF -! -! write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)' - do ig=1,ngrid - 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(:,nlay)-pphi(:,nlay-1))/RG - do l=1,nlay - zlay(:,l)=pphi(:,l)/RG - enddo - do l=1,nlay - deltaz(:,l)=zlev(:,l+1)-zlev(:,l) - enddo - -!----------------------------------------------------------------------- -! Calcul des densites et masses -!----------------------------------------------------------------------- - - 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 - 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 - -!===================================================================== -! Old version of thermcell_plume in thermcell_plume_6A.F90 -! It includes both thermcell_plume_6A and thermcell_plume_5B corresponding -! to the 5B and 6A versions used for CMIP5 and CMIP6. -! The latest was previously named thermcellV1_plume. -! The new thermcell_plume is a clean version (removing obsolete -! options) of thermcell_plume_6A. -! The 3 versions are controled by -! flag_thermals_ed <= 9 thermcell_plume_6A -! <= 19 thermcell_plume_5B -! else thermcell_plume (default 20 for convergence with 6A) -! Fredho -!===================================================================== - - if (iflag_thermals_ed<=9) then -! print*,'THERM NOUVELLE/NOUVELLE Arnaud' - CALL thermcell_plume_6A(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<=19) then -! print*,'THERM RIO et al 2010, version d Arnaud' - CALL thermcell_plume_5B(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) - else - 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) - endif - - if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out - - call test_ltherm(ngrid,nlay,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ') - call test_ltherm(ngrid,nlay,pplay,lmix ,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,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ') - call test_ltherm(ngrid,nlay,pplay,lmin ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin ') - call test_ltherm(ngrid,nlay,pplay,lmix ,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix ') - call test_ltherm(ngrid,nlay,pplay,lmax ,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,pplay,lmin,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry lmin ') -call test_ltherm(ngrid,nlay,pplay,lalim,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_physic (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,pplay,lalim,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ') - call test_ltherm(ngrid,nlay,pplay,lmax ,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) - -#ifdef ISO - ! C Risi: on utilise directement la même routine - do ixt=1,ntraciso - do ll=1,nlay - DO ig=1,ngrid - xtpo_tmp(ig,ll)=xtpo(ixt,ig,ll) - xtzo_tmp(ig,ll)=xtzo(ixt,ig,ll) - enddo - enddo - call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse, & - & xtpo_tmp,xtpdoadj_tmp,xtzo_tmp,lev_out) - do ll=1,nlay - DO ig=1,ngrid - xtpdoadj(ixt,ig,ll)=xtpdoadj_tmp(ig,ll) - enddo - enddo - enddo !do ixt=1,ntraciso -#endif - -#ifdef ISO -#ifdef ISOVERIF - DO ll=1,nlay - DO ig=1,ngrid - if (iso_eau.gt.0) then - call iso_verif_egalite(xtpo(iso_eau,ig,ll), & - & po(ig,ll),'thermcell_main 594') - call iso_verif_egalite(xtpdoadj(iso_eau,ig,ll), & - & pdoadj(ig,ll),'thermcell_main 596') - endif - if (iso_HDO.gt.0) then - call iso_verif_aberrant_encadre(xtpo(iso_hdo,ig,ll) & - & /po(ig,ll),'thermcell_main 610') - endif - enddo - enddo !DO ll=1,nlay - write(*,*) 'thermcell_main 600 tmp: apres thermcell_dq' -#endif -#endif - - - -!------------------------------------------------------------------ -! Calcul de la fraction de l'ascendance -!------------------------------------------------------------------ - do ig=1,ngrid - fraca(ig,1)=0. - fraca(ig,nlay+1)=0. - enddo - do l=2,nlay - do ig=1,ngrid - 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_physic (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 - -!calcul du ratqscdiff - if (prt_level.ge.1) print*,'14e OK convect8' - var=0. - vardiff=0. - ratqsdiff(:,:)=0. - - do l=1,nlay - 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,nlay - 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.) - enddo - enddo - endif - - if (prt_level.ge.1) print*,'thermcell_main FIN OK' - - return - end subroutine thermcell_main - -!============================================================================= -!///////////////////////////////////////////////////////////////////////////// -!============================================================================= - subroutine test_ltherm(ngrid,nlay,pplay,long,ztv,po,ztva, & ! in - & zqla,f_star,zw2,comment) ! in -!============================================================================= - USE thermcell_ini_mod, ONLY: prt_level - IMPLICIT NONE - - integer i, k, ngrid,nlay - real, intent(in), dimension(ngrid,nlay) :: pplay,ztv,po,ztva,zqla - real, intent(in), dimension(ngrid,nlay) :: f_star,zw2 - integer, intent(in), dimension(ngrid) :: long - real seuil - character*21 comment - seuil=0.25 - - if (prt_level.ge.1) THEN - print*,'WARNING !!! TEST ',comment - endif - return - -! test sur la hauteur des thermiques ... - do i=1,ngrid -!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,nlay - 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, & ! in - & therm_tke_max) ! out - USE thermcell_ini_mod, ONLY: prt_level - implicit none - -!======================================================================= -! -! 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 - - real, intent(in) :: ptimestep - real, intent(in), dimension(ngrid,nlay+1) :: fm0,pplev - real, intent(in), dimension(ngrid,nlay) :: entr0 - real, intent(in) :: rg - real, intent(out), dimension(ngrid,nlay) :: therm_tke_max - - real detr0(ngrid,nlay) - real masse0(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) - integer ig,k - - - 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. - - - 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 -