! ! $Id: calltherm.F90 5390 2024-12-05 16:09:25Z ymeurdesoif $ ! !$gpum horizontal klon ngrid MODULE calltherm_mod CONTAINS subroutine calltherm(itap, 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_old=>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 INTEGER, INTENT(IN) :: itap !IM 140508 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),PARAMETER :: modname='calltherm' character (len=80) :: abort_message integer i,k,isplit 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 u_seri(:,:)=u_seri_(:,:) v_seri(:,:)=v_seri_(:,:) t_seri(:,:)=t_seri_(:,:) q_seri(:,:)=q_seri_(:,:) ! 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_old(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*,'<<<<calcul flux ascendant conservatif' fmc_therm=0. do k=1,klev do i=1,klon if (entr_therm(i,k).gt.0.) then fmc_therm(i,k+1)=fmc_therm(i,k)+entr_therm(i,k) else fmc_therm(i,k+1)=fmc_therm(i,k) endif detrc_therm(i,k)=(fmc_therm(i,k+1)-fm_therm(i,k+1)) & & -(fmc_therm(i,k)-fm_therm(i,k)) enddo enddo !**************************************************************** ! calcul de l'humidite dans l'ascendance ! print*,'<<<<calcul de lhumidite dans thermique' !CR:on ne le calcule que pour le cas sec if (iflag_thermals.le.11) then do i=1,klon zqasc(i,1)=q_seri(i,1) do k=2,klev if (fmc_therm(i,k+1).gt.1.e-6) then zqasc(i,k)=(fmc_therm(i,k)*zqasc(i,k-1) & & +entr_therm(i,k)*q_seri(i,k))/fmc_therm(i,k+1) !CR:test on asseche le thermique ! zqasc(i,k)=zqasc(i,k)/2. ! else ! zqasc(i,k)=q_seri(i,k) endif enddo enddo ! calcul de l'eau condensee dans l'ascendance ! print*,'<<<<calcul de leau condensee dans thermique' do i=1,klon do k=1,klev clwcon0(i,k)=zqasc(i,k)-zqsat(i,k) if (clwcon0(i,k).lt.0. .or. & & (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then clwcon0(i,k)=0. endif enddo enddo else do i=1,klon do k=1,klev clwcon0(i,k)=zqla(i,k) if (clwcon0(i,k).lt.0. .or. & & (fm_therm(i,k+1)+detrc_therm(i,k)).lt.1.e-6) then clwcon0(i,k)=0. endif enddo enddo endif !******************************************************************* !jyg Protection contre les temperatures nulles do i=1,klon do k=1,klev if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0. enddo enddo return end subroutine calltherm END MODULE calltherm_mod