Index: LMDZ6/trunk/libf/phylmd/calltherm_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/calltherm_mod.F90 (revision 5513) +++ LMDZ6/trunk/libf/phylmd/calltherm_mod.F90 (revision 5513) @@ -0,0 +1,544 @@ +! +! $Id: calltherm.F90 5390 2024-12-05 16:09:25Z ymeurdesoif $ +! +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*,'<<<