! $Id: calltherm.F90 5158 2024-08-02 12:12:03Z fairhead $ 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 lmdz_print_control, 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 USE lmdz_abort_physic, ONLY: abort_physic USE lmdz_clesphys #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 IMPLICIT NONE 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, 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>=10) THEN PRINT*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion' endif ! tests sur les valeurs negatives de l'eau logexpr0=prt_level>=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)>=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>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',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==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==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==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==11) THEN abort_message = 'cas non prevu dans calltherm' CALL abort_physic (modname,abort_message,1) ELSE IF (iflag_thermals==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>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)>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>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>=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)>=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>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)<50..or.t_seri(i,k)>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>0) PRINT*,'Number of points with q_seri(i,k)<=0 ',nbptspb enddo ! isplit !*************************************************************** ! calcul du flux ascencant conservatif ! PRINT*,'<<<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*,'<<<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*,'<<<