Index: LMDZ6/trunk/libf/phylmd/calltherm.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/calltherm.F90 (revision 5514) +++ (revision ) @@ -1,550 +1,0 @@ -! -! $Id$ -! - subroutine calltherm(dtime & - & ,pplay,paprs,pphi,weak_inversion & - & ,u_seri_,v_seri_,t_seri_,q_seri_,t_env,q_env,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,strig,zcong,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 - USE lmdz_thermcell_alp, ONLY: thermcell_alp - USE lmdz_thermcell_main, ONLY: thermcell_main - USE lmdz_thermcell_old, ONLY: thermcell, thermcell_2002, thermcell_eau, calcul_sec, thermcell_sec -#ifdef ISO - use infotrac_phy, ONLY: ntiso -#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 - USE clesphys_mod_h - USE thermcell_old_mod_h, ONLY: r_aspect_thermals, l_mix_thermals, w2di_thermals - implicit none - - -!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, DIMENSION(klon,klev), INTENT(IN) :: u_seri_,v_seri_ - REAL, DIMENSION(klon,klev), INTENT(IN) :: t_seri_,q_seri_ - REAL, DIMENSION(klon,klev), INTENT(IN) :: t_env,q_env - REAL, DIMENSION(klon,klev) :: u_seri,v_seri - REAL, DIMENSION(klon,klev) :: t_seri,q_seri - REAL, DIMENSION(klon,klev) :: qmemoire - 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 ztv(klon,klev),ztva(klon,klev) - real zpspsk(klon,klev) - real ztla(klon,klev) - real zthl(klon,klev) - real wmax_sec(klon) - real zcong(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),strig(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(ntiso,klon,klev),xtmemoire(ntiso,klon,klev) - REAL d_xt_ajs(ntiso,klon,klev) - real d_xt_the(ntiso,klon,klev) -#ifdef DIAGISO - real q_the(klon,klev) - real xt_the(ntiso,klon,klev) -#endif - real qprec(klon,klev) - integer ixt -#endif - - -!$OMP THREADPRIVATE(first) -!******************************************************** - if (first) then - itap=0 - first=.false. - endif - - u_seri(:,:)=u_seri_(:,:) - v_seri(:,:)=v_seri_(:,:) - t_seri(:,:)=t_seri_(:,:) - q_seri(:,:)=q_seri_(:,:) - -! 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,ntiso - 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_physic('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',ntiso,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,t_env,q_env & - & ,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,zcong & -#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 - & ,zcong,ale,alp,lalim_conv,wght_th & ! out - & ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &! out - & ,n2,s2,strig,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,ntiso - 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',ntiso,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,ntiso - 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',ntiso,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*,'<<<