MODULE cv30_routines_mod !------------------------------------------------------------ ! Parameters for convectL, iflag_con=30: ! (includes - microphysical parameters, ! - parameters that control the rate of approach ! to quasi-equilibrium) ! - noff & minorig (previously in input of convect1) !------------------------------------------------------------ IMPLICIT NONE; PRIVATE PUBLIC sigd, spfac, pbcrit, ptcrit, omtrain, dtovsh, dpbase, dttrig, dtcrit, & tau, beta, alpha, delta, betad, noff, minorig, nl, nlp, nlm, & cv30_param, cv30_prelim, cv30_feed, cv30_undilute1, cv30_trigger, & cv30_compress, cv30_undilute2, cv30_closure, cv30_mixing, cv30_unsat, & cv30_yield, cv30_tracer, cv30_uncompress, cv30_epmax_fn_cape INTEGER noff, minorig, nl, nlp, nlm REAL sigd, spfac REAL pbcrit, ptcrit REAL omtrain REAL dtovsh, dpbase, dttrig REAL dtcrit, tau, beta, alpha REAL delta REAL betad !$OMP THREADPRIVATE(sigd, spfac, pbcrit, ptcrit, omtrain, dtovsh, dpbase, dttrig, dtcrit, & !$OMP tau, beta, alpha, delta, betad, noff, minorig, nl, nlp, nlm) CONTAINS SUBROUTINE cv30_param(nd, delt) USE conema3_mod_h IMPLICIT NONE ! ------------------------------------------------------------ ! Set parameters for convectL for iflag_con = 3 ! ------------------------------------------------------------ ! *** PBCRIT IS THE CRITICAL CLOUD DEPTH (MB) BENEATH WHICH THE *** ! *** PRECIPITATION EFFICIENCY IS ASSUMED TO BE ZERO *** ! *** PTCRIT IS THE CLOUD DEPTH (MB) ABOVE WHICH THE PRECIP. *** ! *** EFFICIENCY IS ASSUMED TO BE UNITY *** ! *** SIGD IS THE FRACTIONAL AREA COVERED BY UNSATURATED DNDRAFT *** ! *** SPFAC IS THE FRACTION OF PRECIPITATION FALLING OUTSIDE *** ! *** OF CLOUD *** ! [TAU: CHARACTERISTIC TIMESCALE USED TO COMPUTE ALPHA & BETA] ! *** ALPHA AND BETA ARE PARAMETERS THAT CONTROL THE RATE OF *** ! *** APPROACH TO QUASI-EQUILIBRIUM *** ! *** (THEIR STANDARD VALUES ARE 1.0 AND 0.96, RESPECTIVELY) *** ! *** (BETA MUST BE LESS THAN OR EQUAL TO 1) *** ! *** DTCRIT IS THE CRITICAL BUOYANCY (K) USED TO ADJUST THE *** ! *** APPROACH TO QUASI-EQUILIBRIUM *** ! *** IT MUST BE LESS THAN 0 *** INTEGER nd REAL delt ! timestep (seconds) ! noff: integer limit for convection (nd-noff) ! minorig: First level of convection ! -- limit levels for convection: noff = 1 minorig = 1 nl = nd - noff nlp = nl + 1 nlm = nl - 1 ! -- "microphysical" parameters: sigd = 0.01 spfac = 0.15 pbcrit = 150.0 ptcrit = 500.0 ! IM cf. FH epmax = 0.993 omtrain = 45.0 ! used also for snow (no disctinction rain/snow) ! -- misc: dtovsh = -0.2 ! dT for overshoot dpbase = -40. ! definition cloud base (400m above LCL) dttrig = 5. ! (loose) condition for triggering ! -- rate of approach to quasi-equilibrium: dtcrit = -2.0 tau = 8000. beta = 1.0 - delt/tau alpha = 1.5E-3*delt/tau ! increase alpha to compensate W decrease: alpha = alpha*1.5 ! -- interface cloud parameterization: delta = 0.01 ! cld ! -- interface with boundary-layer (gust factor): (sb) betad = 10.0 ! original value (from convect 4.3) END SUBROUTINE cv30_param SUBROUTINE cv30_prelim(len, nd, ndp1, t, q, p, ph, lv, cpn, tv, gz, h, hm, & th) USE cvthermo_mod_h IMPLICIT NONE ! ===================================================================== ! --- CALCULATE ARRAYS OF GEOPOTENTIAL, HEAT CAPACITY & STATIC ENERGY ! "ori": from convect4.3 (vectorized) ! "convect3": to be exactly consistent with convect3 ! ===================================================================== ! inputs: INTEGER len, nd, ndp1 REAL t(len, nd), q(len, nd), p(len, nd), ph(len, ndp1) ! outputs: REAL lv(len, nd), cpn(len, nd), tv(len, nd) REAL gz(len, nd), h(len, nd), hm(len, nd) REAL th(len, nd) ! local variables: INTEGER k, i REAL rdcp REAL tvx, tvy ! convect3 REAL cpx(len, nd) ! ori do 110 k=1,nlp DO k = 1, nl ! convect3 DO i = 1, len ! debug lv(i,k)= lv0-clmcpv*(t(i,k)-t0) lv(i, k) = lv0 - clmcpv*(t(i,k)-273.15) cpn(i, k) = cpd*(1.0-q(i,k)) + cpv*q(i, k) cpx(i, k) = cpd*(1.0-q(i,k)) + cl*q(i, k) ! ori tv(i,k)=t(i,k)*(1.0+q(i,k)*epsim1) tv(i, k) = t(i, k)*(1.0+q(i,k)/eps-q(i,k)) rdcp = (rrd*(1.-q(i,k))+q(i,k)*rrv)/cpn(i, k) th(i, k) = t(i, k)*(1000.0/p(i,k))**rdcp END DO END DO ! gz = phi at the full levels (same as p). DO i = 1, len gz(i, 1) = 0.0 END DO ! ori do 140 k=2,nlp DO k = 2, nl ! convect3 DO i = 1, len tvx = t(i, k)*(1.+q(i,k)/eps-q(i,k)) !convect3 tvy = t(i, k-1)*(1.+q(i,k-1)/eps-q(i,k-1)) !convect3 gz(i, k) = gz(i, k-1) + 0.5*rrd*(tvx+tvy) & !convect3 *(p(i,k-1)-p(i,k))/ph(i, k) !convect3 ! ori gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k)) ! ori & *(p(i,k-1)-p(i,k))/ph(i,k) END DO END DO ! h = phi + cpT (dry static energy). ! hm = phi + cp(T-Tbase)+Lq ! ori do 170 k=1,nlp DO k = 1, nl ! convect3 DO i = 1, len h(i, k) = gz(i, k) + cpn(i, k)*t(i, k) hm(i, k) = gz(i, k) + cpx(i, k)*(t(i,k)-t(i,1)) + lv(i, k)*q(i, k) END DO END DO END SUBROUTINE cv30_prelim SUBROUTINE cv30_feed(len, nd, t, q, qs, p, ph, hm, gz, nk, icb, icbmax, & iflag, tnk, qnk, gznk, plcl & #ifdef ISO ,xt,xtnk & #endif ) #ifdef ISO USE infotrac_phy, ONLY: ntraciso=>ntiso #endif IMPLICIT NONE ! ================================================================ ! Purpose: CONVECTIVE FEED ! Main differences with cv_feed: ! - ph added in input ! - here, nk(i)=minorig ! - icb defined differently (plcl compared with ph instead of p) ! Main differences with convect3: ! - we do not compute dplcldt and dplcldr of CLIFT anymore ! - values iflag different (but tests identical) ! - A,B explicitely defined (!...) ! ================================================================ ! inputs: INTEGER len, nd REAL t(len, nd), q(len, nd), qs(len, nd), p(len, nd) REAL hm(len, nd), gz(len, nd) REAL ph(len, nd+1) #ifdef ISO REAL xt(ntraciso,len,nd) #endif ! outputs: INTEGER iflag(len), nk(len), icb(len), icbmax REAL tnk(len), qnk(len), gznk(len), plcl(len) #ifdef ISO REAL xtnk(ntraciso,len) #endif ! local variables: INTEGER i, k #ifdef ISO INTEGER ixt #endif INTEGER ihmin(len) REAL work(len) REAL pnk(len), qsnk(len), rh(len), chi(len) REAL a, b ! convect3 ! ym plcl = 0.0 ! @ !------------------------------------------------------------------- ! @ ! --- Find level of minimum moist static energy ! @ ! --- If level of minimum moist static energy coincides with ! @ ! --- or is lower than minimum allowable parcel origin level, ! @ ! --- set iflag to 6. ! @ !------------------------------------------------------------------- ! @ ! @ do 180 i=1,len ! @ work(i)=1.0e12 ! @ ihmin(i)=nl ! @ 180 continue ! @ do 200 k=2,nlp ! @ do 190 i=1,len ! @ if((hm(i,k).lt.work(i)).AND. ! @ & (hm(i,k).lt.hm(i,k-1)))THEN ! @ work(i)=hm(i,k) ! @ ihmin(i)=k ! @ endif ! @ 190 continue ! @ 200 continue ! @ do 210 i=1,len ! @ ihmin(i)=min(ihmin(i),nlm) ! @ IF(ihmin(i).le.minorig)THEN ! @ iflag(i)=6 ! @ endif ! @ 210 continue ! @ c ! @ !------------------------------------------------------------------- ! @ ! --- Find that model level below the level of minimum moist static ! @ ! --- energy that has the maximum value of moist static energy ! @ !------------------------------------------------------------------- ! @ ! @ do 220 i=1,len ! @ work(i)=hm(i,minorig) ! @ nk(i)=minorig ! @ 220 continue ! @ do 240 k=minorig+1,nl ! @ do 230 i=1,len ! @ if((hm(i,k).gt.work(i)).AND.(k.le.ihmin(i)))THEN ! @ work(i)=hm(i,k) ! @ nk(i)=k ! @ endif ! @ 230 continue ! @ 240 continue ! ------------------------------------------------------------------- ! --- Origin level of ascending parcels for convect3: ! ------------------------------------------------------------------- DO i = 1, len nk(i) = minorig END DO ! ------------------------------------------------------------------- ! --- Check whether parcel level temperature and specific humidity ! --- are reasonable ! ------------------------------------------------------------------- DO i = 1, len IF (((t(i,nk(i))<250.0) .OR. (q(i,nk(i))<=0.0)) & ! @ & .OR.( ! p(i,ihmin(i)).lt.400.0 ! ) ) .AND. (iflag(i)==0)) iflag(i) = 7 END DO ! ------------------------------------------------------------------- ! --- Calculate lifted condensation level of air at parcel origin level ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980) ! ------------------------------------------------------------------- a = 1669.0 ! convect3 b = 122.0 ! convect3 DO i = 1, len IF (iflag(i)/=7) THEN ! modif sb Jun7th 2002 tnk(i) = t(i, nk(i)) qnk(i) = q(i, nk(i)) gznk(i) = gz(i, nk(i)) pnk(i) = p(i, nk(i)) qsnk(i) = qs(i, nk(i)) #ifdef ISO DO ixt=1,ntraciso xtnk(ixt,i) = xt(ixt,i, nk(i)) enddo #endif rh(i) = qnk(i)/qsnk(i) ! ori rh(i)=min(1.0,rh(i)) ! removed for convect3 ! ori chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i)) chi(i) = tnk(i)/(a-b*rh(i)-tnk(i)) ! convect3 plcl(i) = pnk(i)*(rh(i)**chi(i)) IF (((plcl(i)<200.0) .OR. (plcl(i)>=2000.0)) .AND. (iflag(i)==0)) iflag & (i) = 8 END IF ! iflag=7 END DO ! ------------------------------------------------------------------- ! --- Calculate first level above lcl (=icb) ! ------------------------------------------------------------------- ! @ do 270 i=1,len ! @ icb(i)=nlm ! @ 270 continue ! @c ! @ do 290 k=minorig,nl ! @ do 280 i=1,len ! @ if((k.ge.(nk(i)+1)).AND.(p(i,k).lt.plcl(i))) ! @ & icb(i)=min(icb(i),k) ! @ 280 continue ! @ 290 continue ! @c ! @ do 300 i=1,len ! @ if((icb(i).ge.nlm).AND.(iflag(i).EQ.0))iflag(i)=9 ! @ 300 continue DO i = 1, len icb(i) = nlm END DO ! la modification consiste a comparer plcl a ph et non a p: ! icb est defini par : ph(icb)<plcl<ph(icb-1) ! @ do 290 k=minorig,nl DO k = 3, nl - 1 ! modif pour que icb soit sup/egal a 2 DO i = 1, len IF (ph(i,k)<plcl(i)) icb(i) = min(icb(i), k) END DO END DO DO i = 1, len ! @ if((icb(i).ge.nlm).AND.(iflag(i).EQ.0))iflag(i)=9 IF ((icb(i)==nlm) .AND. (iflag(i)==0)) iflag(i) = 9 END DO DO i = 1, len icb(i) = icb(i) - 1 ! icb sup ou egal a 2 END DO ! Compute icbmax. icbmax = 2 DO i = 1, len ! icbmax=max(icbmax,icb(i)) IF (iflag(i)<7) icbmax = max(icbmax, icb(i)) ! sb Jun7th02 END DO END SUBROUTINE cv30_feed SUBROUTINE cv30_undilute1(len, nd, t, q, qs, gz, plcl, p, nk, icb, tp, tvp, & clw, icbs & #ifdef ISO ,xt,xtclw & #endif ) #ifdef ISO USE infotrac_phy, ONLY: ntraciso=>ntiso USE isotopes_mod, ONLY: pxtmelt,pxtice,pxtmin,pxtmax,cond_temp_env, & iso_eau,iso_HDO, ridicule USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall #ifdef ISOTRAC USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall_trac #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: iso_verif_traceur #endif #endif #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: errmax,errmaxrel,Tmin_verif, & iso_verif_egalite_choix, iso_verif_noNaN,iso_verif_aberrant, & iso_verif_egalite,iso_verif_egalite_choix_nostop,iso_verif_positif_nostop, & iso_verif_egalite_nostop,iso_verif_aberrant_nostop,deltaD,iso_verif_noNaN_nostop, & iso_verif_positif,iso_verif_egalite_vect2D #endif #endif USE cvthermo_mod_h IMPLICIT NONE ! ---------------------------------------------------------------- ! Equivalent de TLIFT entre NK et ICB+1 inclus ! Differences with convect4: ! - specify plcl in input ! - icbs is the first level above LCL (may differ from icb) ! - in the iterations, used x(icbs) instead x(icb) ! - many minor differences in the iterations ! - tvp is computed in only one time ! - icbs: first level above Plcl (IMIN de TLIFT) in output ! - if icbs=icb, compute also tp(icb+1),tvp(icb+1) & clw(icb+1) ! ---------------------------------------------------------------- ! inputs: INTEGER len, nd INTEGER nk(len), icb(len) REAL t(len, nd), q(len, nd), qs(len, nd), gz(len, nd) REAL p(len, nd) REAL plcl(len) ! convect3 #ifdef ISO REAL xt(ntraciso,len,nd) #endif ! outputs: REAL tp(len, nd), tvp(len, nd), clw(len, nd) #ifdef ISO REAL xtclw(ntraciso,len,nd) REAL tg_save(len,nd) #endif ! local variables: INTEGER i, k INTEGER icb1(len), icbs(len), icbsmax2 ! convect3 REAL tg, qg, alv, s, ahg, tc, denom, es, rg REAL ah0(len), cpp(len) REAL tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len) REAL qsicb(len) ! convect3 REAL cpinv(len) ! convect3 #ifdef ISO INTEGER ixt REAL zfice(len),zxtliq(ntraciso,len),zxtice(ntraciso,len) REAL q_k(len),clw_k(len),tg_k(len),xt_k(ntraciso,len) !#ifdef ISOVERIF ! integer iso_verif_positif_nostop !#endif #endif ! ------------------------------------------------------------------- ! --- Calculates the lifted parcel virtual temperature at nk, ! --- the actual temperature, and the adiabatic ! --- liquid water content. The procedure is to solve the equation. ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. ! ------------------------------------------------------------------- #ifdef ISOVERIF WRITE(*,*) 'cv30_routine undilute 1 413: entree' #endif DO i = 1, len tnk(i) = t(i, nk(i)) qnk(i) = q(i, nk(i)) gznk(i) = gz(i, nk(i)) ! ori ticb(i)=t(i,icb(i)) ! ori gzicb(i)=gz(i,icb(i)) END DO ! *** Calculate certain parcel quantities, including static energy *** DO i = 1, len ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) + qnk(i)*(lv0-clmcpv*(tnk(i)- & 273.15)) + gznk(i) cpp(i) = cpd*(1.-qnk(i)) + qnk(i)*cpv cpinv(i) = 1./cpp(i) END DO ! *** Calculate lifted parcel quantities below cloud base *** DO i = 1, len !convect3 icb1(i) = min(max(icb(i), 2), nl) ! if icb is below LCL, start loop at ICB+1: ! (icbs est le premier niveau au-dessus du LCL) icbs(i) = icb1(i) !convect3 IF (plcl(i)<p(i,icb1(i))) THEN icbs(i) = min(icbs(i)+1, nl) !convect3 END IF END DO !convect3 DO i = 1, len !convect3 ticb(i) = t(i, icbs(i)) !convect3 gzicb(i) = gz(i, icbs(i)) !convect3 qsicb(i) = qs(i, icbs(i)) !convect3 END DO !convect3 ! Re-compute icbsmax (icbsmax2): !convect3 !convect3 icbsmax2 = 2 !convect3 DO i = 1, len !convect3 icbsmax2 = max(icbsmax2, icbs(i)) !convect3 END DO !convect3 ! initialization outputs: DO k = 1, icbsmax2 ! convect3 DO i = 1, len ! convect3 tp(i, k) = 0.0 ! convect3 tvp(i, k) = 0.0 ! convect3 clw(i, k) = 0.0 ! convect3 #ifdef ISO DO ixt=1,ntraciso xtclw(ixt,i,k) = 0.0 enddo #endif END DO ! convect3 END DO ! convect3 ! tp and tvp below cloud base: DO k = minorig, icbsmax2 - 1 DO i = 1, len tp(i, k) = tnk(i) - (gz(i,k)-gznk(i))*cpinv(i) tvp(i, k) = tp(i, k)*(1.+qnk(i)/eps-qnk(i)) !whole thing (convect3) END DO END DO ! *** Find lifted parcel quantities above cloud base *** DO i = 1, len tg = ticb(i) ! ori qg=qs(i,icb(i)) qg = qsicb(i) ! convect3 ! debug alv=lv0-clmcpv*(ticb(i)-t0) alv = lv0 - clmcpv*(ticb(i)-273.15) ! First iteration. ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3 +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3 s = 1./s ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3 tg = tg + s*(ah0(i)-ahg) ! ori tg=max(tg,35.0) ! debug tc=tg-t0 tc = tg - 273.15 denom = 243.5 + tc denom = max(denom, 1.0) ! convect3 ! ori IF(tc.ge.0.0)THEN es = 6.112*exp(17.67*tc/denom) ! ori else ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) ! ori endif ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) qg = eps*es/(p(i,icbs(i))-es*(1.-eps)) ! qg=max(0.0,qg) ! C Risi ! Second iteration. ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! ori s=1./s ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3 tg = tg + s*(ah0(i)-ahg) ! ori tg=max(tg,35.0) ! debug tc=tg-t0 tc = tg - 273.15 denom = 243.5 + tc denom = max(denom, 1.0) ! convect3 ! ori IF(tc.ge.0.0)THEN es = 6.112*exp(17.67*tc/denom) ! ori else ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) ! ori end if ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) qg = eps*es/(p(i,icbs(i))-es*(1.-eps)) ! qg=max(0.0,qg) ! C Risi alv = lv0 - clmcpv*(ticb(i)-273.15) ! ori c approximation here: ! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i) ! ori & -gz(i,icb(i))-alv*qg)/cpd ! convect3: no approximation: tp(i, icbs(i)) = (ah0(i)-gz(i,icbs(i))-alv*qg)/(cpd+(cl-cpd)*qnk(i)) ! ori clw(i,icb(i))=qnk(i)-qg ! ori clw(i,icb(i))=max(0.0,clw(i,icb(i))) clw(i, icbs(i)) = qnk(i) - qg clw(i, icbs(i)) = max(0.0, clw(i,icbs(i))) rg = qg/(1.-qnk(i)) ! ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi) ! convect3: (qg utilise au lieu du vrai mixing ratio rg) tvp(i, icbs(i)) = tp(i, icbs(i))*(1.+qg/eps-qnk(i)) !whole thing END DO #ifdef ISO ! calcul de zfice DO i=1,len zfice(i) = 1.0-(t(i,icbs(i))-pxtice)/(pxtmelt-pxtice) zfice(i) = MIN(MAX(zfice(i),0.0),1.0) enddo ! calcul de la composition du condensat glace et liquide DO i=1,len clw_k(i)=clw(i,icbs(i)) tg_k(i)=t(i,icbs(i)) DO ixt=1,ntraciso xt_k(ixt,i)=xt(ixt,i,nk(i)) enddo enddo #ifdef ISOVERIF WRITE(*,*) 'cv30_routine undilute1 573: avant condiso' WRITE(*,*) 't(1,1)=',t(1,1) DO i=1,len CALL iso_verif_positif(t(i,icbs(i))-Tmin_verif, & 'cv30_routines 654') enddo IF (iso_HDO.gt.0) THEN DO i=1,len IF (qnk(i).gt.ridicule) THEN CALL iso_verif_aberrant(xt_k(iso_hdo,i)/qnk(i), & 'cv30_routines 576') endif !if (qnk(i).gt.ridicule) THEN enddo endif !if (iso_HDO.gt.0) THEN ! WRITE(*,*) 'i=1, clw_k,qnk=',clw_k(1),qnk(1) #endif CALL condiso_liq_ice_vectall(xt_k(1,1),qnk(1), & clw_k(1),tg_k(1), & zfice(1),zxtice(1,1),zxtliq(1,1),len) #ifdef ISOTRAC #ifdef ISOVERIF WRITE(*,*) 'cv30_routines 658: CALL condiso_liq_ice_vectall_trac' #endif CALL condiso_liq_ice_vectall_trac(xt_k(1,1),qnk(1), & clw_k(1),tg_k(1), & zfice(1),zxtice(1,1),zxtliq(1,1),len) #endif DO i=1,len DO ixt = 1, ntraciso xtclw(ixt,i,icbs(i))= zxtice(ixt,i)+zxtliq(ixt,i) xtclw(ixt,i,icbs(i))=max(0.0,xtclw(ixt,i,icbs(i))) enddo !do ixt=1,niso enddo !do i=1,len #ifdef ISOVERIF WRITE(*,*) 'cv30_routine undilute 1 598: apres condiso' IF (iso_eau.gt.0) THEN DO i=1,len CALL iso_verif_egalite_choix(xtclw(iso_eau,i,icbs(i)), & clw(i,icbs(i)),'cv30_routines 577',errmax,errmaxrel) enddo !do i=1,len endif !if (iso_eau.gt.0) THEN #ifdef ISOTRAC DO i=1,len CALL iso_verif_traceur(xtclw(1,i,k),'cv30_routines 603') enddo #endif #endif #endif ! ori do 380 k=minorig,icbsmax2 ! ori do 370 i=1,len ! ori tvp(i,k)=tvp(i,k)-tp(i,k)*qnk(i) ! ori 370 continue ! ori 380 continue ! -- The following is only for convect3: ! * icbs is the first level above the LCL: ! if plcl<p(icb), then icbs=icb+1 ! if plcl>p(icb), then icbs=icb ! * the routine above computes tvp from minorig to icbs (included). ! * to compute buoybase (in cv3_trigger.F), both tvp(icb) and tvp(icb+1) ! must be known. This is the case if icbs=icb+1, but not if icbs=icb. ! * therefore, in the case icbs=icb, we compute tvp at level icb+1 ! (tvp at other levels will be computed in cv3_undilute2.F) DO i = 1, len ticb(i) = t(i, icb(i)+1) gzicb(i) = gz(i, icb(i)+1) qsicb(i) = qs(i, icb(i)+1) END DO DO i = 1, len tg = ticb(i) qg = qsicb(i) ! convect3 ! debug alv=lv0-clmcpv*(ticb(i)-t0) alv = lv0 - clmcpv*(ticb(i)-273.15) ! First iteration. ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3 +alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! convect3 s = 1./s ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3 tg = tg + s*(ah0(i)-ahg) ! ori tg=max(tg,35.0) ! debug tc=tg-t0 tc = tg - 273.15 denom = 243.5 + tc denom = max(denom, 1.0) ! convect3 ! ori IF(tc.ge.0.0)THEN es = 6.112*exp(17.67*tc/denom) ! ori else ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) ! ori endif ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps)) ! qg=max(0.0,qg) ! C Risi ! Second iteration. ! ori s=cpd+alv*alv*qg/(rrv*ticb(i)*ticb(i)) ! ori s=1./s ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*ticb(i)+alv*qg+gzicb(i) ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gzicb(i) ! convect3 tg = tg + s*(ah0(i)-ahg) ! ori tg=max(tg,35.0) ! debug tc=tg-t0 tc = tg - 273.15 denom = 243.5 + tc denom = max(denom, 1.0) ! convect3 ! ori IF(tc.ge.0.0)THEN es = 6.112*exp(17.67*tc/denom) ! ori else ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) ! ori end if ! ori qg=eps*es/(p(i,icb(i))-es*(1.-eps)) qg = eps*es/(p(i,icb(i)+1)-es*(1.-eps)) ! qg=max(0.0,qg) ! C Risi alv = lv0 - clmcpv*(ticb(i)-273.15) ! ori c approximation here: ! ori tp(i,icb(i))=(ah0(i)-(cl-cpd)*qnk(i)*ticb(i) ! ori & -gz(i,icb(i))-alv*qg)/cpd ! convect3: no approximation: tp(i, icb(i)+1) = (ah0(i)-gz(i,icb(i)+1)-alv*qg)/(cpd+(cl-cpd)*qnk(i)) ! ori clw(i,icb(i))=qnk(i)-qg ! ori clw(i,icb(i))=max(0.0,clw(i,icb(i))) clw(i, icb(i)+1) = qnk(i) - qg clw(i, icb(i)+1) = max(0.0, clw(i,icb(i)+1)) rg = qg/(1.-qnk(i)) ! ori tvp(i,icb(i))=tp(i,icb(i))*(1.+rg*epsi) ! convect3: (qg utilise au lieu du vrai mixing ratio rg) tvp(i, icb(i)+1) = tp(i, icb(i)+1)*(1.+qg/eps-qnk(i)) !whole thing END DO #ifdef ISO DO i=1,len zfice(i) = 1.0-(t(i,icb(i)+1)-pxtice)/(pxtmelt-pxtice) zfice(i) = MIN(MAX(zfice(i),0.0),1.0) ! CALL calcul_zfice(tp(i,icb(i)+1),zfice) enddo !do i=1,len DO i=1,len clw_k(i)=clw(i,icb(i)+1) tg_k(i)=t(i,icb(i)+1) #ifdef ISOVERIF CALL iso_verif_positif(tg_k(i)-Tmin_verif,'cv30_routines 750') #endif DO ixt=1,ntraciso xt_k(ixt,i)=xt(ixt,i,nk(i)) enddo enddo !do i=1,len #ifdef ISOVERIF WRITE(*,*) 'cv30_routines 739: avant condiso' IF (iso_HDO.gt.0) THEN DO i=1,len CALL iso_verif_aberrant(xt_k(iso_hdo,i)/qnk(i), & 'cv30_routines 725') enddo endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC DO i=1,len CALL iso_verif_traceur(xtclw(1,i,k),'cv30_routines 738') enddo #endif #endif CALL condiso_liq_ice_vectall(xt_k(1,1),qnk(1), & clw_k(1),tg_k(1), & zfice(1),zxtice(1,1),zxtliq(1,1),len) #ifdef ISOTRAC CALL condiso_liq_ice_vectall_trac(xt_k(1,1),qnk(1), & clw_k(1),tg_k(1), & zfice(1),zxtice(1,1),zxtliq(1,1),len) #endif DO i=1,len DO ixt = 1, ntraciso xtclw(ixt,i,icb(i)+1)=zxtice(ixt,i)+zxtliq(ixt,i) xtclw(ixt,i,icb(i)+1)=max(0.0,xtclw(ixt,i,icb(i)+1)) enddo !do ixt = 1, niso enddo !do i=1,len #ifdef ISOVERIF !WRITE(*,*) 'DEBUG ISO B' DO i=1,len IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xtclw(iso_eau,i,icb(i)+1), & clw(i,icb(i)+1),'cv30_routines 708',errmax,errmaxrel) endif ! if (iso_eau.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur(xtclw(1,i,icb(i)+1), & 'cv30_routines 760') #endif enddo !do i=1,len !WRITE(*,*) 'FIN DEBUG ISO B' #endif #endif END SUBROUTINE cv30_undilute1 SUBROUTINE cv30_trigger(len, nd, icb, plcl, p, th, tv, tvp, pbase, buoybase, & iflag, sig, w0) IMPLICIT NONE ! ------------------------------------------------------------------- ! --- TRIGGERING ! - computes the cloud base ! - triggering (crude in this version) ! - relaxation of sig and w0 when no convection ! Caution1: if no convection, we set iflag=4 ! (it used to be 0 in convect3) ! Caution2: at this stage, tvp (and thus buoy) are know up ! through icb only! ! -> the buoyancy below cloud base not (yet) set to the cloud base buoyancy ! ------------------------------------------------------------------- ! input: INTEGER len, nd INTEGER icb(len) REAL plcl(len), p(len, nd) REAL th(len, nd), tv(len, nd), tvp(len, nd) ! output: REAL pbase(len), buoybase(len) ! input AND output: INTEGER iflag(len) REAL sig(len, nd), w0(len, nd) ! local variables: INTEGER i, k REAL tvpbase, tvbase, tdif, ath, ath1 ! *** set cloud base buoyancy at (plcl+dpbase) level buoyancy DO i = 1, len pbase(i) = plcl(i) + dpbase tvpbase = tvp(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ & (p(i,icb(i))-p(i,icb(i)+1)) + tvp(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/( & p(i,icb(i))-p(i,icb(i)+1)) tvbase = tv(i, icb(i))*(pbase(i)-p(i,icb(i)+1))/ & (p(i,icb(i))-p(i,icb(i)+1)) + tv(i, icb(i)+1)*(p(i,icb(i))-pbase(i))/(p & (i,icb(i))-p(i,icb(i)+1)) buoybase(i) = tvpbase - tvbase END DO ! *** make sure that column is dry adiabatic between the surface *** ! *** and cloud base, and that lifted air is positively buoyant *** ! *** at cloud base *** ! *** if not, return to calling program after resetting *** ! *** sig(i) and w0(i) *** ! oct3 do 200 i=1,len ! oct3 ! oct3 tdif = buoybase(i) ! oct3 ath1 = th(i,1) ! oct3 ath = th(i,icb(i)-1) - dttrig ! oct3 ! oct3 if (tdif.lt.dtcrit .OR. ath.gt.ath1) THEN ! oct3 do 60 k=1,nl ! oct3 sig(i,k) = beta*sig(i,k) - 2.*alpha*tdif*tdif ! oct3 sig(i,k) = AMAX1(sig(i,k),0.0) ! oct3 w0(i,k) = beta*w0(i,k) ! oct3 60 continue ! oct3 iflag(i)=4 ! pour version vectorisee ! oct3c convect3 iflag(i)=0 ! oct3cccc RETURN ! oct3 endif ! oct3 ! oct3200 continue ! -- oct3: on reecrit la boucle 200 (pour la vectorisation) DO k = 1, nl DO i = 1, len tdif = buoybase(i) ath1 = th(i, 1) ath = th(i, icb(i)-1) - dttrig IF (tdif<dtcrit .OR. ath>ath1) THEN sig(i, k) = beta*sig(i, k) - 2.*alpha*tdif*tdif sig(i, k) = amax1(sig(i,k), 0.0) w0(i, k) = beta*w0(i, k) iflag(i) = 4 ! pour version vectorisee ! convect3 iflag(i)=0 END IF END DO END DO ! fin oct3 -- END SUBROUTINE cv30_trigger SUBROUTINE cv30_compress(len, nloc, ncum, nd, ntra, iflag1, nk1, icb1, icbs1, & plcl1, tnk1, qnk1, gznk1, pbase1, buoybase1, t1, q1, qs1, u1, v1, gz1, & th1, tra1, h1, lv1, cpn1, p1, ph1, tv1, tp1, tvp1, clw1, sig1, w01, & iflag, nk, icb, icbs, plcl, tnk, qnk, gznk, pbase, buoybase, t, q, qs, u, & v, gz, th, tra, h, lv, cpn, p, ph, tv, tp, tvp, clw, sig, w0 & #ifdef ISO ,xtnk1,xt1,xtclw1 & ,xtnk,xt,xtclw & #endif ) USE print_control_mod, ONLY: lunout #ifdef ISO USE infotrac_phy, ONLY: ntraciso=>ntiso USE isotopes_mod, ONLY: essai_convergence, iso_eau,iso_HDO #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: errmax,errmaxrel,Tmin_verif, & iso_verif_egalite_choix, iso_verif_noNaN,iso_verif_aberrant, & iso_verif_egalite,iso_verif_egalite_choix_nostop,iso_verif_positif_nostop, & iso_verif_egalite_nostop,iso_verif_aberrant_nostop,deltaD,iso_verif_noNaN_nostop, & iso_verif_positif,iso_verif_egalite_vect2D #endif #endif IMPLICIT NONE ! inputs: INTEGER len, ncum, nd, ntra, nloc INTEGER iflag1(len), nk1(len), icb1(len), icbs1(len) REAL plcl1(len), tnk1(len), qnk1(len), gznk1(len) REAL pbase1(len), buoybase1(len) REAL t1(len, nd), q1(len, nd), qs1(len, nd), u1(len, nd), v1(len, nd) REAL gz1(len, nd), h1(len, nd), lv1(len, nd), cpn1(len, nd) REAL p1(len, nd), ph1(len, nd+1), tv1(len, nd), tp1(len, nd) REAL tvp1(len, nd), clw1(len, nd) REAL th1(len, nd) REAL sig1(len, nd), w01(len, nd) REAL tra1(len, nd, ntra) #ifdef ISO !integer niso REAL xt1(ntraciso,len,nd), xtclw1(ntraciso,len,nd) REAL xtnk1(ntraciso,len) #endif ! outputs: ! en fait, on a nloc=len pour l'instant (cf cv_driver) INTEGER iflag(nloc), nk(nloc), icb(nloc), icbs(nloc) REAL plcl(nloc), tnk(nloc), qnk(nloc), gznk(nloc) REAL pbase(nloc), buoybase(nloc) REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), u(nloc, nd), v(nloc, nd) REAL gz(nloc, nd), h(nloc, nd), lv(nloc, nd), cpn(nloc, nd) REAL p(nloc, nd), ph(nloc, nd+1), tv(nloc, nd), tp(nloc, nd) REAL tvp(nloc, nd), clw(nloc, nd) REAL th(nloc, nd) REAL sig(nloc, nd), w0(nloc, nd) REAL tra(nloc, nd, ntra) #ifdef ISO REAL xt(ntraciso,nloc,nd), xtclw(ntraciso,nloc,nd) REAL xtnk(ntraciso,nloc) #endif ! local variables: INTEGER i, k, nn, j #ifdef ISO INTEGER ixt #endif CHARACTER (LEN=20) :: modname = 'cv30_compress' CHARACTER (LEN=80) :: abort_message #ifdef ISO ! initialisation des champs compresses: DO k=1,nd DO i=1,nloc IF (essai_convergence) THEN else q(i,k)=0.0 clw(i,k)=0.0 ! mise en commentaire le 5 avril pour verif ! convergence endif !f (negation(essai_convergence)) THEN DO ixt=1,ntraciso xt(ixt,i,k)=0.0 xtclw(ixt,i,k)=0.0 enddo !do ixt=1,niso enddo !do i=1,len enddo !do k=1,nd ! WRITE(*,*) 'q(1,1),xt(iso_eau,1,1)=',q(1,1),xt(iso_eau,1,1) #endif DO k = 1, nl + 1 nn = 0 DO i = 1, len IF (iflag1(i)==0) THEN nn = nn + 1 sig(nn, k) = sig1(i, k) w0(nn, k) = w01(i, k) t(nn, k) = t1(i, k) q(nn, k) = q1(i, k) qs(nn, k) = qs1(i, k) u(nn, k) = u1(i, k) v(nn, k) = v1(i, k) gz(nn, k) = gz1(i, k) h(nn, k) = h1(i, k) lv(nn, k) = lv1(i, k) cpn(nn, k) = cpn1(i, k) p(nn, k) = p1(i, k) ph(nn, k) = ph1(i, k) tv(nn, k) = tv1(i, k) tp(nn, k) = tp1(i, k) tvp(nn, k) = tvp1(i, k) clw(nn, k) = clw1(i, k) th(nn, k) = th1(i, k) #ifdef ISO DO ixt = 1, ntraciso xt(ixt,nn,k)=xt1(ixt,i,k) xtclw(ixt,nn,k)=xtclw1(ixt,i,k) enddo ! WRITE(*,*) 'nn,i,k,q(nn,k),xt(iso_eau,nn,k)=', & ! & nn,i,k,q(nn, k),xt(ixt,nn,k) #endif END IF END DO END DO ! do 121 j=1,ntra ! do 111 k=1,nd ! nn=0 ! do 101 i=1,len ! IF(iflag1(i).EQ.0)THEN ! nn=nn+1 ! tra(nn,k,j)=tra1(i,k,j) ! END IF ! 101 continue ! 111 continue ! 121 continue IF (nn/=ncum) THEN WRITE (lunout, *) 'strange! nn not equal to ncum: ', nn, ncum abort_message = '' CALL abort_physic(modname, abort_message, 1) END IF nn = 0 DO i = 1, len IF (iflag1(i)==0) THEN nn = nn + 1 pbase(nn) = pbase1(i) buoybase(nn) = buoybase1(i) plcl(nn) = plcl1(i) tnk(nn) = tnk1(i) qnk(nn) = qnk1(i) gznk(nn) = gznk1(i) nk(nn) = nk1(i) icb(nn) = icb1(i) icbs(nn) = icbs1(i) iflag(nn) = iflag1(i) #ifdef ISO DO ixt=1,ntraciso xtnk(ixt,nn) = xtnk1(ixt,i) enddo #endif END IF END DO #ifdef ISO #ifdef ISOVERIF IF (iso_eau.gt.0) THEN DO k = 1, nd DO i = 1, nloc !WRITE(*,*) 'i,k=',i,k CALL iso_verif_egalite_choix(xtclw(iso_eau,i,k),clw(i,k), & 'compress 973',errmax,errmaxrel) CALL iso_verif_egalite_choix(xt(iso_eau,i,k),q(i,k), & 'compress 975',errmax,errmaxrel) enddo enddo endif !if (iso_eau.gt.0) THEN DO k = 1, nd DO i = 1, nn CALL iso_verif_positif(q(i,k),'compress 1004') enddo enddo #endif #endif END SUBROUTINE cv30_compress SUBROUTINE cv30_undilute2(nloc, ncum, nd, icb, icbs, nk, tnk, qnk, gznk, t, & q, qs, gz, p, h, tv, lv, pbase, buoybase, plcl, inb, tp, tvp, clw, hp, & ep, sigp, buoy & #ifdef ISO ,xtnk,xt,xtclw & #endif ) ! epmax_cape: ajout arguments USE conema3_mod_h #ifdef ISO USE infotrac_phy, ONLY: ntraciso=>ntiso USE isotopes_mod, ONLY: pxtmelt,pxtice,pxtmin,pxtmax,cond_temp_env, iso_eau,iso_HDO USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall #ifdef ISOTRAC USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall_trac #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: iso_verif_traceur #endif #endif #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: errmax,errmaxrel,Tmin_verif,Tmax_verif, & iso_verif_egalite_choix, iso_verif_noNaN,iso_verif_aberrant, & iso_verif_egalite,iso_verif_egalite_choix_nostop,iso_verif_positif_nostop, & iso_verif_egalite_nostop,iso_verif_aberrant_nostop,deltaD,iso_verif_noNaN_nostop, & iso_verif_positif,iso_verif_egalite_vect2D #endif #endif USE cvthermo_mod_h IMPLICIT NONE ! --------------------------------------------------------------------- ! Purpose: ! FIND THE REST OF THE LIFTED PARCEL TEMPERATURES ! & ! COMPUTE THE PRECIPITATION EFFICIENCIES AND THE ! FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD ! & ! FIND THE LEVEL OF NEUTRAL BUOYANCY ! Main differences convect3/convect4: ! - icbs (input) is the first level above LCL (may differ from icb) ! - many minor differences in the iterations ! - condensed water not removed from tvp in convect3 ! - vertical profile of buoyancy computed here (use of buoybase) ! - the determination of inb is different ! - no inb1, ONLY inb in output ! --------------------------------------------------------------------- ! inputs: INTEGER ncum, nd, nloc INTEGER icb(nloc), icbs(nloc), nk(nloc) REAL t(nloc, nd), q(nloc, nd), qs(nloc, nd), gz(nloc, nd) REAL p(nloc, nd) REAL tnk(nloc), qnk(nloc), gznk(nloc) REAL lv(nloc, nd), tv(nloc, nd), h(nloc, nd) REAL pbase(nloc), buoybase(nloc), plcl(nloc) ! outputs: INTEGER inb(nloc) REAL tp(nloc, nd), tvp(nloc, nd), clw(nloc, nd) REAL ep(nloc, nd), sigp(nloc, nd), hp(nloc, nd) REAL buoy(nloc, nd) ! local variables: INTEGER i, k REAL tg, qg, ahg, alv, s, tc, es, denom, rg, tca, elacrit REAL by, defrac, pden REAL ah0(nloc), cape(nloc), capem(nloc), byp(nloc) LOGICAL lcape(nloc) #ifdef ISO REAL xt(ntraciso,nloc,nd), xtclw(ntraciso,nloc,nd) REAL xtnk(ntraciso,nloc) ! real xtep(ntraciso,nloc,nd) ! le 7 mai: on supprime xtep, car pas besoin ! la chute de precip ne fractionne pas. INTEGER ixt REAL zfice(nloc),zxtliq(ntraciso,nloc),zxtice(ntraciso,nloc) REAL clw_k(nloc),tg_k(nloc) #ifdef ISOVERIF REAL qg_save(nloc,nd) ! inout !integer iso_verif_positif_nostop #endif #endif ! ===================================================================== ! --- SOME INITIALIZATIONS ! ===================================================================== DO k = 1, nl DO i = 1, ncum ep(i, k) = 0.0 sigp(i, k) = spfac clw(i,k)=0.0 ! C Risi END DO END DO ! ===================================================================== ! --- FIND THE REST OF THE LIFTED PARCEL TEMPERATURES ! ===================================================================== ! --- The procedure is to solve the equation. ! cp*tp+L*qp+phi=cp*tnk+L*qnk+gznk. ! *** Calculate certain parcel quantities, including static energy *** DO i = 1, ncum ah0(i) = (cpd*(1.-qnk(i))+cl*qnk(i))*tnk(i) & ! debug & ! +qnk(i)*(lv0-clmcpv*(tnk(i)-t0))+gznk(i) +qnk(i)*(lv0-clmcpv*(tnk(i)-273.15)) + gznk(i) END DO ! *** Find lifted parcel quantities above cloud base *** DO k = minorig + 1, nl DO i = 1, ncum ! ori IF(k.ge.(icb(i)+1))THEN IF (k>=(icbs(i)+1)) THEN ! convect3 tg = t(i, k) qg = qs(i, k) ! debug alv=lv0-clmcpv*(t(i,k)-t0) alv = lv0 - clmcpv*(t(i,k)-273.15) ! First iteration. ! ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k)) s = cpd*(1.-qnk(i)) + cl*qnk(i) & ! convect3 +alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! convect3 s = 1./s ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k) ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3 tg = tg + s*(ah0(i)-ahg) ! ori tg=max(tg,35.0) ! debug tc=tg-t0 tc = tg - 273.15 denom = 243.5 + tc denom = max(denom, 1.0) ! convect3 ! ori IF(tc.ge.0.0)THEN es = 6.112*exp(17.67*tc/denom) ! ori else ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) ! ori endif qg = eps*es/(p(i,k)-es*(1.-eps)) ! qg=max(0.0,qg) ! C Risi ! Second iteration. ! ori s=cpd+alv*alv*qg/(rrv*t(i,k)*t(i,k)) ! ori s=1./s ! ori ahg=cpd*tg+(cl-cpd)*qnk(i)*t(i,k)+alv*qg+gz(i,k) ahg = cpd*tg + (cl-cpd)*qnk(i)*tg + alv*qg + gz(i, k) ! convect3 tg = tg + s*(ah0(i)-ahg) ! ori tg=max(tg,35.0) ! debug tc=tg-t0 tc = tg - 273.15 denom = 243.5 + tc denom = max(denom, 1.0) ! convect3 ! ori IF(tc.ge.0.0)THEN es = 6.112*exp(17.67*tc/denom) ! ori else ! ori es=exp(23.33086-6111.72784/tg+0.15215*log(tg)) ! ori endif qg = eps*es/(p(i,k)-es*(1.-eps)) ! qg=max(0.0,qg) ! C Risi ! debug alv=lv0-clmcpv*(t(i,k)-t0) alv = lv0 - clmcpv*(t(i,k)-273.15) ! PRINT*,'cpd dans convect2 ',cpd ! PRINT*,'tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd' ! PRINT*,tp(i,k),ah0(i),cl,cpd,qnk(i),t(i,k),gz(i,k),alv,qg,cpd ! ori c approximation here: ! ori ! tp(i,k)=(ah0(i)-(cl-cpd)*qnk(i)*t(i,k)-gz(i,k)-alv*qg)/cpd ! convect3: no approximation: tp(i, k) = (ah0(i)-gz(i,k)-alv*qg)/(cpd+(cl-cpd)*qnk(i)) clw(i, k) = qnk(i) - qg clw(i, k) = max(0.0, clw(i,k)) rg = qg/(1.-qnk(i)) ! ori tvp(i,k)=tp(i,k)*(1.+rg*epsi) ! convect3: (qg utilise au lieu du vrai mixing ratio rg): tvp(i, k) = tp(i, k)*(1.+qg/eps-qnk(i)) ! whole thing END IF END DO #ifdef ISO ! calcul de zfice DO i=1,ncum zfice(i) = 1.0-(t(i,k)-pxtice)/(pxtmelt-pxtice) zfice(i) = MIN(MAX(zfice(i),0.0),1.0) enddo DO i=1,ncum clw_k(i)=clw(i,k) tg_k(i)=t(i,k) enddo !do i=1,ncum #ifdef ISOVERIF !WRITE(*,*) 'cv30_routine 1259: avant condiso' IF (iso_HDO.gt.0) THEN DO i=1,ncum CALL iso_verif_aberrant(xtnk(iso_hdo,i)/qnk(i), & 'cv30_routines 1231') enddo endif !if (iso_HDO.gt.0) THEN IF (iso_eau.gt.0) THEN DO i=1,ncum CALL iso_verif_egalite(xtnk(iso_eau,i),qnk(i), & 'cv30_routines 1373') enddo endif !if (iso_HDO.gt.0) THEN DO i=1,ncum IF ((iso_verif_positif_nostop(qnk(i)-clw_k(i), & 'cv30_routines 1275').EQ.1).OR. & (iso_verif_positif_nostop(tg_k(i)-Tmin_verif, & 'cv30_routines 1297a').EQ.1).OR. & (iso_verif_positif_nostop(Tmax_verif-tg_k(i), & 'cv30_routines 1297b').EQ.1)) THEN WRITE(*,*) 'i,k,qnk,clw_k=',i,k,qnk(i),clw_k(i) WRITE(*,*) 'tg,t,qg=',tg_k(i),t(i,k),qg_save(i,k) WRITE(*,*) 'icbs(i)=',icbs(i) stop endif ! if ((iso_verif_positif_nostop enddo !do i=1,ncum #ifdef ISOTRAC DO i=1,ncum CALL iso_verif_traceur(xtnk(1,i),'cv30_routines 1251') enddo !do i=1,ncum #endif #endif CALL condiso_liq_ice_vectall(xtnk(1,1),qnk(1), & clw_k(1),tg_k(1), & zfice(1),zxtice(1,1),zxtliq(1,1),ncum) #ifdef ISOTRAC #ifdef ISOVERIF WRITE(*,*) 'cv30_routines 1283: condiso pour traceurs' #endif CALL condiso_liq_ice_vectall_trac(xtnk(1,1),qnk(1), & clw_k(1),tg_k(1), & zfice(1),zxtice(1,1),zxtliq(1,1),ncum) #endif DO i=1,ncum DO ixt=1,ntraciso xtclw(ixt,i,k)=zxtice(ixt,i)+zxtliq(ixt,i) xtclw(ixt,i,k)=max(0.0,xtclw(ixt,i,k)) enddo !do ixt=1,niso enddo !do i=1,ncum #ifdef ISOVERIF IF (iso_eau.gt.0) THEN DO i=1,ncum CALL iso_verif_egalite_choix(xtclw(iso_eau,i,k), & clw(i,k),'cv30_routines 1223',errmax,errmaxrel) enddo endif !if (iso_eau.gt.0) THEN #ifdef ISOTRAC DO i=1,ncum CALL iso_verif_traceur(xtclw(1,i,k),'cv30_routines 1275') enddo #endif #endif #endif END DO ! ===================================================================== ! --- SET THE PRECIPITATION EFFICIENCIES AND THE FRACTION OF ! --- PRECIPITATION FALLING OUTSIDE OF CLOUD ! --- THESE MAY BE FUNCTIONS OF TP(I), P(I) AND CLW(I) ! ===================================================================== ! ori do 320 k=minorig+1,nl DO k = 1, nl ! convect3 DO i = 1, ncum pden = ptcrit - pbcrit ep(i, k) = (plcl(i)-p(i,k)-pbcrit)/pden*epmax ep(i, k) = amax1(ep(i,k), 0.0) ep(i, k) = amin1(ep(i,k), epmax) sigp(i, k) = spfac ! ori IF(k.ge.(nk(i)+1))THEN ! ori tca=tp(i,k)-t0 ! ori IF(tca.ge.0.0)THEN ! ori elacrit=elcrit ! ori else ! ori elacrit=elcrit*(1.0-tca/tlcrit) ! ori endif ! ori elacrit=max(elacrit,0.0) ! ori ep(i,k)=1.0-elacrit/max(clw(i,k),1.0e-8) ! ori ep(i,k)=max(ep(i,k),0.0 ) ! ori ep(i,k)=min(ep(i,k),1.0 ) ! ori sigp(i,k)=sigs ! ori endif END DO END DO ! ===================================================================== ! --- CALCULATE VIRTUAL TEMPERATURE AND LIFTED PARCEL ! --- VIRTUAL TEMPERATURE ! ===================================================================== ! dans convect3, tvp est calcule en une seule fois, et sans retirer ! l'eau condensee (~> reversible CAPE) ! ori do 340 k=minorig+1,nl ! ori do 330 i=1,ncum ! ori IF(k.ge.(icb(i)+1))THEN ! ori tvp(i,k)=tvp(i,k)*(1.0-qnk(i)+ep(i,k)*clw(i,k)) ! oric PRINT*,'i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k)' ! oric PRINT*, i,k,tvp(i,k),qnk(i),ep(i,k),clw(i,k) ! ori endif ! ori 330 continue ! ori 340 continue ! ori do 350 i=1,ncum ! ori tvp(i,nlp)=tvp(i,nl)-(gz(i,nlp)-gz(i,nl))/cpd ! ori 350 continue DO i = 1, ncum ! convect3 tp(i, nlp) = tp(i, nl) ! convect3 END DO ! convect3 ! ===================================================================== ! --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only): ! ===================================================================== ! -- this is for convect3 only: ! first estimate of buoyancy: DO i = 1, ncum DO k = 1, nl buoy(i, k) = tvp(i, k) - tv(i, k) END DO END DO ! set buoyancy=buoybase for all levels below base ! for safety, set buoy(icb)=buoybase DO i = 1, ncum DO k = 1, nl IF ((k>=icb(i)) .AND. (k<=nl) .AND. (p(i,k)>=pbase(i))) THEN buoy(i, k) = buoybase(i) END IF END DO ! IM cf. CRio/JYG 270807 buoy(icb(i),k)=buoybase(i) buoy(i, icb(i)) = buoybase(i) END DO ! -- end convect3 ! ===================================================================== ! --- FIND THE FIRST MODEL LEVEL (INB) ABOVE THE PARCEL'S ! --- LEVEL OF NEUTRAL BUOYANCY ! ===================================================================== ! -- this is for convect3 only: DO i = 1, ncum inb(i) = nl - 1 END DO DO i = 1, ncum DO k = 1, nl - 1 IF ((k>=icb(i)) .AND. (buoy(i,k)<dtovsh)) THEN inb(i) = min(inb(i), k) END IF END DO END DO ! -- end convect3 ! ori do 510 i=1,ncum ! ori cape(i)=0.0 ! ori capem(i)=0.0 ! ori inb(i)=icb(i)+1 ! ori inb1(i)=inb(i) ! ori 510 continue ! Originial Code ! do 530 k=minorig+1,nl-1 ! do 520 i=1,ncum ! IF(k.ge.(icb(i)+1))THEN ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) ! byp=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) ! cape(i)=cape(i)+by ! IF(by.ge.0.0)inb1(i)=k+1 ! IF(cape(i).gt.0.0)THEN ! inb(i)=k+1 ! capem(i)=cape(i) ! END IF ! END IF ! 520 continue ! 530 continue ! do 540 i=1,ncum ! byp=(tvp(i,nl)-tv(i,nl))*dph(i,nl)/p(i,nl) ! cape(i)=capem(i)+byp ! defrac=capem(i)-cape(i) ! defrac=max(defrac,0.001) ! frac(i)=-cape(i)/defrac ! frac(i)=min(frac(i),1.0) ! frac(i)=max(frac(i),0.0) ! 540 continue ! K Emanuel fix ! CALL zilch(byp,ncum) ! do 530 k=minorig+1,nl-1 ! do 520 i=1,ncum ! IF(k.ge.(icb(i)+1))THEN ! by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) ! cape(i)=cape(i)+by ! IF(by.ge.0.0)inb1(i)=k+1 ! IF(cape(i).gt.0.0)THEN ! inb(i)=k+1 ! capem(i)=cape(i) ! byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) ! END IF ! END IF ! 520 continue ! 530 continue ! do 540 i=1,ncum ! inb(i)=max(inb(i),inb1(i)) ! cape(i)=capem(i)+byp(i) ! defrac=capem(i)-cape(i) ! defrac=max(defrac,0.001) ! frac(i)=-cape(i)/defrac ! frac(i)=min(frac(i),1.0) ! frac(i)=max(frac(i),0.0) ! 540 continue ! J Teixeira fix ! ori CALL zilch(byp,ncum) ! ori do 515 i=1,ncum ! ori lcape(i)=.TRUE. ! ori 515 continue ! ori do 530 k=minorig+1,nl-1 ! ori do 520 i=1,ncum ! ori IF(cape(i).lt.0.0)lcape(i)=.FALSE. ! ori if((k.ge.(icb(i)+1)).AND.lcape(i))THEN ! ori by=(tvp(i,k)-tv(i,k))*dph(i,k)/p(i,k) ! ori byp(i)=(tvp(i,k+1)-tv(i,k+1))*dph(i,k+1)/p(i,k+1) ! ori cape(i)=cape(i)+by ! ori IF(by.ge.0.0)inb1(i)=k+1 ! ori IF(cape(i).gt.0.0)THEN ! ori inb(i)=k+1 ! ori capem(i)=cape(i) ! ori endif ! ori endif ! ori 520 continue ! ori 530 continue ! ori do 540 i=1,ncum ! ori cape(i)=capem(i)+byp(i) ! ori defrac=capem(i)-cape(i) ! ori defrac=max(defrac,0.001) ! ori frac(i)=-cape(i)/defrac ! ori frac(i)=min(frac(i),1.0) ! ori frac(i)=max(frac(i),0.0) ! ori 540 continue ! ===================================================================== ! --- CALCULATE LIQUID WATER STATIC ENERGY OF LIFTED PARCEL ! ===================================================================== ! ym do i=1,ncum*nlp ! ym hp(i,1)=h(i,1) ! ym enddo DO k = 1, nlp DO i = 1, ncum hp(i, k) = h(i, k) END DO END DO DO k = minorig + 1, nl DO i = 1, ncum IF ((k>=icb(i)) .AND. (k<=inb(i))) THEN hp(i, k) = h(i, nk(i)) + (lv(i,k)+(cpd-cpv)*t(i,k))*ep(i, k)*clw(i, k & ) END IF END DO END DO END SUBROUTINE cv30_undilute2 SUBROUTINE cv30_closure(nloc, ncum, nd, icb, inb, pbase, p, ph, tv, buoy, & sig, w0, cape, m) USE cvthermo_mod_h IMPLICIT NONE ! =================================================================== ! --- CLOSURE OF CONVECT3 ! vectorization: S. Bony ! =================================================================== ! input: INTEGER ncum, nd, nloc INTEGER icb(nloc), inb(nloc) REAL pbase(nloc) REAL p(nloc, nd), ph(nloc, nd+1) REAL tv(nloc, nd), buoy(nloc, nd) ! input/output: REAL sig(nloc, nd), w0(nloc, nd) ! output: REAL cape(nloc) REAL m(nloc, nd) ! local variables: INTEGER i, j, k, icbmax REAL deltap, fac, w, amu REAL dtmin(nloc, nd), sigold(nloc, nd) ! ------------------------------------------------------- ! -- Initialization ! ------------------------------------------------------- DO k = 1, nl DO i = 1, ncum m(i, k) = 0.0 END DO END DO ! ------------------------------------------------------- ! -- Reset sig(i) and w0(i) for i>inb and i<icb ! ------------------------------------------------------- ! update sig and w0 above LNB: DO k = 1, nl - 1 DO i = 1, ncum IF ((inb(i)<(nl-1)) .AND. (k>=(inb(i)+1))) THEN sig(i, k) = beta*sig(i, k) + 2.*alpha*buoy(i, inb(i))*abs(buoy(i,inb( & i))) sig(i, k) = amax1(sig(i,k), 0.0) w0(i, k) = beta*w0(i, k) END IF END DO END DO ! compute icbmax: icbmax = 2 DO i = 1, ncum icbmax = max(icbmax, icb(i)) END DO ! update sig and w0 below cloud base: DO k = 1, icbmax DO i = 1, ncum IF (k<=icb(i)) THEN sig(i, k) = beta*sig(i, k) - 2.*alpha*buoy(i, icb(i))*buoy(i, icb(i)) sig(i, k) = amax1(sig(i,k), 0.0) w0(i, k) = beta*w0(i, k) END IF END DO END DO ! IF(inb.lt.(nl-1))THEN ! do 85 i=inb+1,nl-1 ! sig(i)=beta*sig(i)+2.*alpha*buoy(inb)* ! 1 abs(buoy(inb)) ! sig(i)=amax1(sig(i),0.0) ! w0(i)=beta*w0(i) ! 85 continue ! end if ! do 87 i=1,icb ! sig(i)=beta*sig(i)-2.*alpha*buoy(icb)*buoy(icb) ! sig(i)=amax1(sig(i),0.0) ! w0(i)=beta*w0(i) ! 87 continue ! ------------------------------------------------------------- ! -- Reset fractional areas of updrafts and w0 at initial time ! -- and after 10 time steps of no convection ! ------------------------------------------------------------- DO k = 1, nl - 1 DO i = 1, ncum IF (sig(i,nd)<1.5 .OR. sig(i,nd)>12.0) THEN sig(i, k) = 0.0 w0(i, k) = 0.0 END IF END DO END DO ! ------------------------------------------------------------- ! -- Calculate convective available potential energy (cape), ! -- vertical velocity (w), fractional area covered by ! -- undilute updraft (sig), and updraft mass flux (m) ! ------------------------------------------------------------- DO i = 1, ncum cape(i) = 0.0 END DO ! compute dtmin (minimum buoyancy between ICB and given level k): DO i = 1, ncum DO k = 1, nl dtmin(i, k) = 100.0 END DO END DO DO i = 1, ncum DO k = 1, nl DO j = minorig, nl IF ((k>=(icb(i)+1)) .AND. (k<=inb(i)) .AND. (j>=icb(i)) .AND. (j<=(k- & 1))) THEN dtmin(i, k) = amin1(dtmin(i,k), buoy(i,j)) END IF END DO END DO END DO ! the interval on which cape is computed starts at pbase : DO k = 1, nl DO i = 1, ncum IF ((k>=(icb(i)+1)) .AND. (k<=inb(i))) THEN deltap = min(pbase(i), ph(i,k-1)) - min(pbase(i), ph(i,k)) cape(i) = cape(i) + rrd*buoy(i, k-1)*deltap/p(i, k-1) cape(i) = amax1(0.0, cape(i)) sigold(i, k) = sig(i, k) ! dtmin(i,k)=100.0 ! do 97 j=icb(i),k-1 ! mauvaise vectorisation ! dtmin(i,k)=AMIN1(dtmin(i,k),buoy(i,j)) ! 97 continue sig(i, k) = beta*sig(i, k) + alpha*dtmin(i, k)*abs(dtmin(i,k)) sig(i, k) = amax1(sig(i,k), 0.0) sig(i, k) = amin1(sig(i,k), 0.01) fac = amin1(((dtcrit-dtmin(i,k))/dtcrit), 1.0) w = (1.-beta)*fac*sqrt(cape(i)) + beta*w0(i, k) amu = 0.5*(sig(i,k)+sigold(i,k))*w m(i, k) = amu*0.007*p(i, k)*(ph(i,k)-ph(i,k+1))/tv(i, k) w0(i, k) = w END IF END DO END DO DO i = 1, ncum w0(i, icb(i)) = 0.5*w0(i, icb(i)+1) m(i, icb(i)) = 0.5*m(i, icb(i)+1)*(ph(i,icb(i))-ph(i,icb(i)+1))/ & (ph(i,icb(i)+1)-ph(i,icb(i)+2)) sig(i, icb(i)) = sig(i, icb(i)+1) sig(i, icb(i)-1) = sig(i, icb(i)) END DO ! cape=0.0 ! do 98 i=icb+1,inb ! deltap = min(pbase,ph(i-1))-min(pbase,ph(i)) ! cape=cape+rrd*buoy(i-1)*deltap/p(i-1) ! dcape=rrd*buoy(i-1)*deltap/p(i-1) ! dlnp=deltap/p(i-1) ! cape=amax1(0.0,cape) ! sigold=sig(i) ! dtmin=100.0 ! do 97 j=icb,i-1 ! dtmin=amin1(dtmin,buoy(j)) ! 97 continue ! sig(i)=beta*sig(i)+alpha*dtmin*abs(dtmin) ! sig(i)=amax1(sig(i),0.0) ! sig(i)=amin1(sig(i),0.01) ! fac=amin1(((dtcrit-dtmin)/dtcrit),1.0) ! w=(1.-beta)*fac*sqrt(cape)+beta*w0(i) ! amu=0.5*(sig(i)+sigold)*w ! m(i)=amu*0.007*p(i)*(ph(i)-ph(i+1))/tv(i) ! w0(i)=w ! 98 continue ! w0(icb)=0.5*w0(icb+1) ! m(icb)=0.5*m(icb+1)*(ph(icb)-ph(icb+1))/(ph(icb+1)-ph(icb+2)) ! sig(icb)=sig(icb+1) ! sig(icb-1)=sig(icb) END SUBROUTINE cv30_closure SUBROUTINE cv30_mixing(nloc, ncum, nd, na, ntra, icb, nk, inb, ph, t, rr, rs, & u, v, tra, h, lv, qnk, hp, tv, tvp, ep, clw, m, sig, ment, qent, uent, & vent, sij, elij, ments, qents, traent & #ifdef ISO ,xt,xtnk,xtclw & ,xtent,xtelij & #endif ) #ifdef ISO USE infotrac_phy, ONLY: ntraciso=>ntiso,niso,itZonIso USE isotopes_mod, ONLY: pxtmelt,pxtice,pxtmin,pxtmax, iso_eau,iso_HDO, & ridicule USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: errmax,errmaxrel,Tmin_verif,deltalim, & iso_verif_egalite_choix,iso_verif_aberrant_choix, iso_verif_noNaN, & iso_verif_aberrant, & iso_verif_egalite,iso_verif_egalite_choix_nostop,iso_verif_positif_nostop, & iso_verif_egalite_nostop,iso_verif_aberrant_nostop,deltaD,iso_verif_noNaN_nostop, & iso_verif_positif,iso_verif_egalite_vect2D #endif #ifdef ISOTRAC USE isotrac_mod, ONLY: option_tmin,option_traceurs,seuil_tag_tmin, & & option_cond,index_zone,izone_cond,index_iso USE isotrac_routines_mod, ONLY: iso_recolorise_condensation USE isotopes_routines_mod, ONLY: condiso_liq_ice_vectall_trac #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: iso_verif_trac17_q_deltad,iso_verif_traceur, & & iso_verif_traceur_justmass #endif #endif #endif USE cvthermo_mod_h IMPLICIT NONE ! --------------------------------------------------------------------- ! a faire: ! - changer rr(il,1) -> qnk(il) ! - vectorisation de la partie normalisation des flux (do 789...) ! --------------------------------------------------------------------- ! inputs: INTEGER ncum, nd, na, ntra, nloc INTEGER icb(nloc), inb(nloc), nk(nloc) REAL sig(nloc, nd) REAL qnk(nloc) REAL ph(nloc, nd+1) REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd) REAL u(nloc, nd), v(nloc, nd) REAL tra(nloc, nd, ntra) ! input of convect3 REAL lv(nloc, na), h(nloc, na), hp(nloc, na) REAL tv(nloc, na), tvp(nloc, na), ep(nloc, na), clw(nloc, na) REAL m(nloc, na) ! input of convect3 #ifdef ISO REAL xt(ntraciso,nloc,na), xtclw(ntraciso,nloc,na) REAL tg_save(nloc,nd) REAL xtnk(ntraciso,nloc) ! real xtep(ntraciso,nloc,na) #endif ! outputs: REAL ment(nloc, na, na), qent(nloc, na, na) REAL uent(nloc, na, na), vent(nloc, na, na) REAL sij(nloc, na, na), elij(nloc, na, na) REAL traent(nloc, nd, nd, ntra) REAL ments(nloc, nd, nd), qents(nloc, nd, nd) REAL sigij(nloc, nd, nd) #ifdef ISO REAL xtent(ntraciso,nloc,nd,nd) REAL xtelij(ntraciso,nloc,nd,nd) #endif ! local variables: INTEGER i, j, k, il, im, jm INTEGER num1, num2 INTEGER nent(nloc, na) REAL rti, bf2, anum, denom, dei, altem, cwat, stemp, qp REAL alt, smid, sjmin, sjmax, delp, delm REAL asij(nloc), smax(nloc), scrit(nloc) REAL asum(nloc, nd), bsum(nloc, nd), csum(nloc, nd) REAL wgh REAL zm(nloc, na) LOGICAL lwork(nloc) #ifdef ISO INTEGER ixt REAL xtrti(ntraciso,nloc) REAL xtres(ntraciso) ! on ajoute la dimension nloc a xtrti pour verifs dans les tags: 5 fev ! 2010 REAL zfice(nloc),zxtliq(ntraciso,nloc),zxtice(ntraciso,nloc) ! real xt_reduit(ntraciso) ! LOGICAL negation !#ifdef ISOVERIF ! integer iso_verif_positif_nostop !#endif #endif ! ===================================================================== ! --- INITIALIZE VARIOUS ARRAYS USED IN THE COMPUTATIONS ! ===================================================================== #ifdef ISO #ifdef ISOVERIF WRITE(*,*) 'cv30_routines 1820: entree dans cv3_mixing' IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_vect2D( & xtclw,clw, & 'cv30_mixing 1841',ntraciso,nloc,na) endif #endif #endif ! ori do 360 i=1,ncum*nlp DO j = 1, nl DO i = 1, ncum nent(i, j) = 0 ! in convect3, m is computed in cv3_closure ! ori m(i,1)=0.0 END DO END DO ! ori do 400 k=1,nlp ! ori do 390 j=1,nlp DO j = 1, nl DO k = 1, nl DO i = 1, ncum qent(i, k, j) = rr(i, j) uent(i, k, j) = u(i, j) vent(i, k, j) = v(i, j) elij(i, k, j) = 0.0 ! ym ment(i,k,j)=0.0 ! ym sij(i,k,j)=0.0 END DO END DO END DO #ifdef ISO DO j=1,nd DO k=1,nd DO i=1,ncum DO ixt =1,ntraciso xtent(ixt,i,k,j)=xt(ixt,i,j) xtelij(ixt,i,k,j)=0.0 enddo !do ixt =1,niso ! on initialise mieux que ca qent et elij, meme si au final les ! valeurs en nd=nl+1 ne sont pas utilisees qent(i,k,j)=rr(i,j) elij(i,k,j)=0.0 enddo !do i=1,ncum enddo !do k=1,nl enddo !do j=1,nl #endif ! ym ment(1:ncum, 1:nd, 1:nd) = 0.0 sij(1:ncum, 1:nd, 1:nd) = 0.0 ! do k=1,ntra ! do j=1,nd ! instead nlp ! do i=1,nd ! instead nlp ! do il=1,ncum ! traent(il,i,j,k)=tra(il,j,k) ! enddo ! enddo ! enddo ! enddo zm(:, :) = 0. ! ===================================================================== ! --- CALCULATE ENTRAINED AIR MASS FLUX (ment), TOTAL WATER MIXING ! --- RATIO (QENT), TOTAL CONDENSED WATER (elij), AND MIXING ! --- FRACTION (sij) ! ===================================================================== DO i = minorig + 1, nl DO j = minorig, nl DO il = 1, ncum IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (j>=(icb(il)- & 1)) .AND. (j<=inb(il))) THEN rti = rr(il, 1) - ep(il, i)*clw(il, i) bf2 = 1. + lv(il, j)*lv(il, j)*rs(il, j)/(rrv*t(il,j)*t(il,j)*cpd) anum = h(il, j) - hp(il, i) + (cpv-cpd)*t(il, j)*(rti-rr(il,j)) denom = h(il, i) - hp(il, i) + (cpd-cpv)*(rr(il,i)-rti)*t(il, j) dei = denom IF (abs(dei)<0.01) dei = 0.01 sij(il, i, j) = anum/dei sij(il, i, i) = 1.0 altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - rs(il, j) altem = altem/bf2 cwat = clw(il, j)*(1.-ep(il,j)) stemp = sij(il, i, j) IF ((stemp<0.0 .OR. stemp>1.0 .OR. altem>cwat) .AND. j>i) THEN anum = anum - lv(il, j)*(rti-rs(il,j)-cwat*bf2) denom = denom + lv(il, j)*(rr(il,i)-rti) IF (abs(denom)<0.01) denom = 0.01 sij(il, i, j) = anum/denom altem = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti - & rs(il, j) altem = altem - (bf2-1.)*cwat END IF IF (sij(il,i,j)>0.0 .AND. sij(il,i,j)<0.95) THEN qent(il, i, j) = sij(il, i, j)*rr(il, i) + (1.-sij(il,i,j))*rti uent(il, i, j) = sij(il, i, j)*u(il, i) + & (1.-sij(il,i,j))*u(il, nk(il)) vent(il, i, j) = sij(il, i, j)*v(il, i) + & (1.-sij(il,i,j))*v(il, nk(il)) ! !!! do k=1,ntra ! !!! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k) ! !!! : +(1.-sij(il,i,j))*tra(il,nk(il),k) ! !!! END DO elij(il, i, j) = altem elij(il, i, j) = amax1(0.0, elij(il,i,j)) ment(il, i, j) = m(il, i)/(1.-sij(il,i,j)) nent(il, i) = nent(il, i) + 1 END IF sij(il, i, j) = amax1(0.0, sij(il,i,j)) sij(il, i, j) = amin1(1.0, sij(il,i,j)) END IF ! new END DO #ifdef ISO #ifdef ISOVERIF !WRITE(*,*) 'cv30_routines tmp 2078' #endif DO il=1,ncum zfice(il) = 1.0-(t(il,j)-pxtice)/(pxtmelt-pxtice) zfice(il) = MIN(MAX(zfice(il),0.0),1.0) IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND. & (j.ge.(icb(il)-1)).AND.(j.le.inb(il)))THEN DO ixt=1,ntraciso ! xtrti(ixt)=xt(ixt,il,1)-xtep(ixt,il,i)*xtclw(ixt,il,i) ! le 7 mai: on supprime xtep xtrti(ixt,il)=xt(ixt,il,1)-ep(il,i)*xtclw(ixt,il,i) enddo IF(sij(il,i,j).gt.0.0.AND.sij(il,i,j).lt.0.95)THEN ! temperature of condensation (within mixtures): ! tcond(il)=t(il,j) ! : + ( sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti ! : - elij(il,i,j) - rs(il,j) ) ! : / ( cpd*(bf2-1.0)/lv(il,j) ) DO ixt = 1, ntraciso ! total mixing ratio in the mixtures before precipitation: xtent(ixt,il,i,j)=sij(il,i,j)*xt(ixt,il,i) & +(1.-sij(il,i,j))*xtrti(ixt,il) enddo !do ixt = 1, ntraciso endif !IF(sij(il,i,j).gt.0.0.AND.sij(il,i,j).lt.0.95)THEN endif !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND. enddo !do il=1,ncum CALL condiso_liq_ice_vectall(xtent(1,1,i,j),qent(1,i,j), & elij(1,i,j), & t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum) #ifdef ISOTRAC CALL condiso_liq_ice_vectall_trac(xtent(1,1,i,j),qent(1,i,j), & elij(1,i,j), & t(1,j),zfice(1),zxtice(1,1),zxtliq(1,1),ncum) #ifdef ISOVERIF DO il=1,ncum CALL iso_verif_traceur(xt(1,il,i),'cv30_routines 1967') IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND. & (j.ge.(icb(il)-1)).AND.(j.le.inb(il)))THEN CALL iso_verif_traceur(xtrti(1,il),'cv30_routines 1968') endif !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND. CALL iso_verif_traceur(xtent(1,il,i,j),'cv30_routines 1969') enddo !do il=1,ncum #endif #endif DO il=1,ncum DO ixt = 1, ntraciso xtelij(ixt,il,i,j)=zxtice(ixt,il)+zxtliq(ixt,il) enddo !do ixt = 1, ntraciso enddo !do il=1,ncum #ifdef ISOVERIF IF ((j.EQ.15).AND.(i.EQ.15)) THEN il=2722 IF (il.le.ncum) THEN WRITE(*,*) 'cv30_routines tmp 2194, il,i,j=',il,i,j WRITE(*,*) 'qent,elij=',qent(il,i,j),elij(il,i,j) WRITE(*,*) 'tgsave,zfice=',t(il,j),zfice(il) WRITE(*,*) 'deltaDqent=',deltaD(xtent(iso_HDO,il,i,j)/qent(il,i,j)) WRITE(*,*) 'deltaDelij=',deltaD(xtelij(iso_HDO,il,i,j)/elij(il,i,j)) WRITE(*,*) 'deltaDice=',deltaD(zxtice(iso_HDO,il)/(zfice(il)*elij(il,i,j))) WRITE(*,*) 'deltaDliq=',deltaD(zxtliq(iso_HDO,il)/(1.0-zfice(il)*elij(il,i,j))) endif endif #endif #ifdef ISOTRAC ! WRITE(*,*) 'cv30_routines tmp 1987,option_traceurs=', ! : option_traceurs IF (option_tmin.ge.1) THEN DO il=1,ncum ! WRITE(*,*) 'cv3 tmp 1991 il,i,j,xtent(:,il,i,j),', ! : 'tcond(il),rs(il,j)=', ! : il,i,j,xtent(:,il,i,j),tcond(il),rs(il,j) ! colorier la vapeur residuelle selon temperature de ! condensation, et le condensat en un tag spEcifique IF ((elij(il,i,j).gt.0.0).AND.(qent(il,i,j).gt.0.0)) THEN IF (option_traceurs.EQ.17) THEN CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), & xtent(1,il,i,j),xtelij(1,il,i,j),t(1,j), & 0.0,xtres, & seuil_tag_tmin) else !if (option_traceurs.EQ.17) THEN ! WRITE(*,*) 'cv3 2002: il,i,j =',il,i,j CALL iso_recolorise_condensation(qent(il,i,j),elij(il,i,j), & xtent(1,il,i,j),xtelij(1,il,i,j),rs(il,j),0.0,xtres, & seuil_tag_tmin) endif !if (option_traceurs.EQ.17) THEN DO ixt=1+niso,ntraciso xtent(ixt,il,i,j)=xtres(ixt) enddo endif !if (cond.gt.0.0) THEN enddo !do il=1,ncum #ifdef ISOVERIF DO il=1,ncum CALL iso_verif_traceur(xtent(1,il,i,j),'cv30_routines 1996') CALL iso_verif_traceur(xtelij(1,il,i,j),'cv30_routines 1997') CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), & 'cv30_routines 2042') enddo !do il=1,ncum #endif endif !if (option_tmin.ge.1) THEN #endif ! fractionation: #ifdef ISOVERIF ! WRITE(*,*) 'cv30_routines 2050: avant condiso' DO il=1,ncum IF ((i.ge.icb(il)).AND.(i.le.inb(il)).AND. & (j.ge.(icb(il)-1)).AND.(j.le.inb(il))) THEN IF (sij(il,i,j).gt.0.0.AND.sij(il,i,j).lt.0.95) THEN IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,j), & qent(il,i,j),'cv30_routines 1889',errmax,errmaxrel) CALL iso_verif_egalite_choix(xtelij(iso_eau,il,i,j), & elij(il,i,j),'cv30_routines 1890',errmax,errmaxrel) endif IF (iso_HDO.gt.0) THEN CALL iso_verif_aberrant_choix(xt(iso_HDO,il,i),rr(il,i), & ridicule,deltalim,'cv30_routines 1997') CALL iso_verif_aberrant_choix( & xtent(iso_HDO,il,i,j),qent(il,i,j), & ridicule,deltalim,'cv30_routines 1931') CALL iso_verif_aberrant_choix( & xtelij(iso_HDO,il,i,j),elij(il,i,j), & ridicule,deltalim,'cv30_routines 1993') endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC ! WRITE(*,*) 'cv30_routines tmp 2039 il=',il CALL iso_verif_traceur(xtent(1,il,i,j), & 'cv30_routines 2031') CALL iso_verif_traceur(xtelij(1,il,i,j), & 'cv30_routines 2033') #endif endif !IF(sij(il,i,j).gt.0.0.AND.sij(il,i,j).lt.0.95)THEN endif !IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND. enddo !do il=1,ncum #endif ! WRITE(*,*) 'cv30_routine tmp 1984: cond=',elij(il,i,j) #endif END DO ! do k=1,ntra ! do j=minorig,nl ! do il=1,ncum ! IF( (i.ge.icb(il)).AND.(i.le.inb(il)).AND. ! : (j.ge.(icb(il)-1)).AND.(j.le.inb(il)))THEN ! traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k) ! : +(1.-sij(il,i,j))*tra(il,nk(il),k) ! END IF ! enddo ! enddo ! enddo ! *** if no air can entrain at level i assume that updraft detrains ! *** ! *** at that level and calculate detrained air flux and properties ! *** ! @ do 170 i=icb(il),inb(il) DO il = 1, ncum IF ((i>=icb(il)) .AND. (i<=inb(il)) .AND. (nent(il,i)==0)) THEN ! @ IF(nent(il,i).EQ.0)THEN ment(il, i, i) = m(il, i) qent(il, i, i) = rr(il, nk(il)) - ep(il, i)*clw(il, i) uent(il, i, i) = u(il, nk(il)) vent(il, i, i) = v(il, nk(il)) elij(il, i, i) = clw(il, i) ! MAF sij(il,i,i)=1.0 sij(il, i, i) = 0.0 #ifdef ISO DO ixt = 1, ntraciso xtent(ixt,il,i,i)=xt(ixt,il,nk(il))-ep(il,i)*xtclw(ixt,il,i) ! xtent(ixt,il,i,i)=xt(ixt,il,nk(il))-xtep(ixt,il,i)*xtclw(ixt,il,i) ! le 7 mai: on supprime xtep xtelij(ixt,il,i,i)=xtclw(ixt,il,i) ! rq: ne sera pas utilise ensuite enddo !do ixt = 1, ntraciso #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xtelij(iso_eau,il,i,i), & elij(il,i,i),'cv30_mixing 2117',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN #endif #ifdef ISOTRAC IF (option_tmin.ge.1) THEN ! colorier la vapeur residuelle selon temperature de ! condensation, et le condensat en un tag specifique ! WRITE(*,*) 'cv3 tmp 2095 il,i,j,xtent(:,il,i,j)=', ! : il,i,j,xtent(:,il,i,j) IF ((elij(il,i,i).gt.0.0).AND.(qent(il,i,i).gt.0.0)) THEN IF (option_traceurs.EQ.17) THEN CALL iso_recolorise_condensation(qent(il,i,i), & elij(il,i,i), & xt(1,il,nk(il)),xtclw(1,il,i),t(il,i),ep(il,i), & xtres, & seuil_tag_tmin) else !if (option_traceurs.EQ.17) THEN CALL iso_recolorise_condensation(qent(il,i,i), & elij(il,i,i), & xt(1,il,nk(il)),xtclw(1,il,i),rs(il,i),ep(il,i), & xtres, & seuil_tag_tmin) endif !if (option_traceurs.EQ.17) THEN DO ixt=1+niso,ntraciso xtent(ixt,il,i,i)=xtres(ixt) enddo #ifdef ISOVERIF DO ixt=1,niso CALL iso_verif_egalite_choix(xtres(ixt),xtent(ixt,il,i,i), & 'cv30_routines 2102',errmax,errmaxrel) CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), & 'cv30_routines 2154') enddo #endif endif !if (cond.gt.0.0) THEN #ifdef ISOVERIF CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,i), & qent(il,i,i),'cv30_routines 2103',errmax,errmaxrel) CALL iso_verif_traceur(xtent(1,il,i,i),'cv30_routines 2095') CALL iso_verif_traceur(xtelij(1,il,i,i),'cv30_routines 2096') #endif endif !if (option_tmin.ge.1) THEN #endif #endif END IF END DO END DO ! do j=1,ntra ! do i=minorig+1,nl ! do il=1,ncum ! if (i.ge.icb(il) .AND. i.le.inb(il) .AND. nent(il,i).EQ.0) THEN ! traent(il,i,i,j)=tra(il,nk(il),j) ! END IF ! enddo ! enddo ! enddo DO j = minorig, nl DO i = minorig, nl DO il = 1, ncum IF ((j>=(icb(il)-1)) .AND. (j<=inb(il)) .AND. (i>=icb(il)) .AND. (i<= & inb(il))) THEN sigij(il, i, j) = sij(il, i, j) END IF END DO END DO END DO ! @ enddo ! @170 continue ! ===================================================================== ! --- NORMALIZE ENTRAINED AIR MASS FLUXES ! --- TO REPRESENT EQUAL PROBABILITIES OF MIXING ! ===================================================================== ! ym CALL zilch(asum,ncum*nd) ! ym CALL zilch(bsum,ncum*nd) ! ym CALL zilch(csum,ncum*nd) CALL zilch(asum, nloc*nd) CALL zilch(csum, nloc*nd) CALL zilch(csum, nloc*nd) DO il = 1, ncum lwork(il) = .FALSE. END DO DO i = minorig + 1, nl num1 = 0 DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il)) num1 = num1 + 1 END DO IF (num1<=0) GO TO 789 DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il)) THEN lwork(il) = (nent(il,i)/=0) qp = rr(il, 1) - ep(il, i)*clw(il, i) anum = h(il, i) - hp(il, i) - lv(il, i)*(qp-rs(il,i)) + & (cpv-cpd)*t(il, i)*(qp-rr(il,i)) denom = h(il, i) - hp(il, i) + lv(il, i)*(rr(il,i)-qp) + & (cpd-cpv)*t(il, i)*(rr(il,i)-qp) IF (abs(denom)<0.01) denom = 0.01 scrit(il) = anum/denom alt = qp - rs(il, i) + scrit(il)*(rr(il,i)-qp) IF (scrit(il)<=0.0 .OR. alt<=0.0) scrit(il) = 1.0 smax(il) = 0.0 asij(il) = 0.0 END IF END DO DO j = nl, minorig, -1 num2 = 0 DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( & il)-1) .AND. j<=inb(il) .AND. lwork(il)) num2 = num2 + 1 END DO IF (num2<=0) GO TO 175 DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. j>=(icb( & il)-1) .AND. j<=inb(il) .AND. lwork(il)) THEN IF (sij(il,i,j)>1.0E-16 .AND. sij(il,i,j)<0.95) THEN wgh = 1.0 IF (j>i) THEN sjmax = amax1(sij(il,i,j+1), smax(il)) sjmax = amin1(sjmax, scrit(il)) smax(il) = amax1(sij(il,i,j), smax(il)) sjmin = amax1(sij(il,i,j-1), smax(il)) sjmin = amin1(sjmin, scrit(il)) IF (sij(il,i,j)<(smax(il)-1.0E-16)) wgh = 0.0 smid = amin1(sij(il,i,j), scrit(il)) ELSE sjmax = amax1(sij(il,i,j+1), scrit(il)) smid = amax1(sij(il,i,j), scrit(il)) sjmin = 0.0 IF (j>1) sjmin = sij(il, i, j-1) sjmin = amax1(sjmin, scrit(il)) END IF delp = abs(sjmax-smid) delm = abs(sjmin-smid) asij(il) = asij(il) + wgh*(delp+delm) ment(il, i, j) = ment(il, i, j)*(delp+delm)*wgh END IF END IF END DO 175 END DO DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN asij(il) = amax1(1.0E-16, asij(il)) asij(il) = 1.0/asij(il) asum(il, i) = 0.0 bsum(il, i) = 0.0 csum(il, i) = 0.0 END IF END DO DO j = minorig, nl DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( & il)-1) .AND. j<=inb(il)) THEN ment(il, i, j) = ment(il, i, j)*asij(il) END IF END DO END DO DO j = minorig, nl DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( & il)-1) .AND. j<=inb(il)) THEN asum(il, i) = asum(il, i) + ment(il, i, j) ment(il, i, j) = ment(il, i, j)*sig(il, j) bsum(il, i) = bsum(il, i) + ment(il, i, j) END IF END DO END DO DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il)) THEN bsum(il, i) = amax1(bsum(il,i), 1.0E-16) bsum(il, i) = 1.0/bsum(il, i) END IF END DO DO j = minorig, nl DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( & il)-1) .AND. j<=inb(il)) THEN ment(il, i, j) = ment(il, i, j)*asum(il, i)*bsum(il, i) END IF END DO END DO DO j = minorig, nl DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. j>=(icb( & il)-1) .AND. j<=inb(il)) THEN csum(il, i) = csum(il, i) + ment(il, i, j) END IF END DO END DO DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. lwork(il) .AND. & csum(il,i)<m(il,i)) THEN nent(il, i) = 0 ment(il, i, i) = m(il, i) qent(il, i, i) = rr(il, 1) - ep(il, i)*clw(il, i) uent(il, i, i) = u(il, nk(il)) vent(il, i, i) = v(il, nk(il)) elij(il, i, i) = clw(il, i) ! MAF sij(il,i,i)=1.0 sij(il, i, i) = 0.0 #ifdef ISO DO ixt = 1, ntraciso ! xtent(ixt,il,i,i)=xt(ixt,il,1)-xtep(ixt,il,i)*xtclw(ixt,il,i) xtent(ixt,il,i,i)=xt(ixt,il,1)-ep(il,i)*xtclw(ixt,il,i) xtelij(ixt,il,i,i)=xtclw(ixt,il,i) ! rq: ne sera pas utilise ensuite enddo #endif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xtelij(iso_eau,il,i,i), & elij(il,i,i),'cv30_mixing 2354',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN #endif #ifdef ISOTRAC IF (option_tmin.ge.1) THEN ! colorier la vapeur residuelle selon temperature de ! condensation, et le condensat en un tag specifique ! WRITE(*,*) 'cv3 tmp 2314 il,i,j,xtent(:,il,i,j)=', ! : il,i,j,xtent(:,il,i,j) IF ((elij(il,i,i).gt.0.0).AND.(qent(il,i,i).gt.0.0)) THEN IF (option_traceurs.EQ.17) THEN CALL iso_recolorise_condensation(qent(il,i,i), & elij(il,i,i), & xt(1,il,1),xtclw(1,il,i),t(il,i),ep(il,i), & xtres, & seuil_tag_tmin) else !if (option_traceurs.EQ.17) THEN CALL iso_recolorise_condensation(qent(il,i,i), & elij(il,i,i), & xt(1,il,1),xtclw(1,il,i),rs(il,i),ep(il,i), & xtres, & seuil_tag_tmin) endif ! if (option_traceurs.EQ.17) THEN DO ixt=1+niso,ntraciso xtent(ixt,il,i,i)=xtres(ixt) enddo #ifdef ISOVERIF DO ixt=1,niso CALL iso_verif_egalite_choix(xtres(ixt),xtent(ixt,il,i,i), & 'cv30_routines 2318',errmax,errmaxrel) CALL iso_verif_trac17_q_deltaD(xtent(1,il,i,j), & 'cv30_routines 2383') enddo #endif endif !if (cond.gt.0.0) THEN #ifdef ISOVERIF CALL iso_verif_egalite_choix(xtent(iso_eau,il,i,i), & qent(il,i,i),'cv30_routines 2321',errmax,errmaxrel) CALL iso_verif_traceur(xtent(1,il,i,i),'cv30_routines 2322') CALL iso_verif_traceur(xtelij(1,il,i,i),'cv30_routines 2323') #endif endif !if (option_tmin.ge.1) THEN #endif END IF END DO ! il ! do j=1,ntra ! do il=1,ncum ! if ( i.ge.icb(il) .AND. i.le.inb(il) .AND. lwork(il) ! : .AND. csum(il,i).lt.m(il,i) ) THEN ! traent(il,i,i,j)=tra(il,nk(il),j) ! END IF ! enddo ! enddo 789 END DO ! MAF: renormalisation de MENT DO jm = 1, nd DO im = 1, nd DO il = 1, ncum zm(il, im) = zm(il, im) + (1.-sij(il,im,jm))*ment(il, im, jm) END DO END DO END DO DO jm = 1, nd DO im = 1, nd DO il = 1, ncum IF (zm(il,im)/=0.) THEN ment(il, im, jm) = ment(il, im, jm)*m(il, im)/zm(il, im) END IF END DO END DO END DO DO jm = 1, nd DO im = 1, nd DO il = 1, ncum qents(il, im, jm) = qent(il, im, jm) ments(il, im, jm) = ment(il, im, jm) END DO END DO END DO #ifdef ISO !c--debug #ifdef ISOVERIF DO im = 1, nd DO jm = 1, nd DO il = 1, ncum IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xtelij(iso_eau,il,im,jm), & elij(il,im,jm),'cv30_mixing 2110',errmax,errmaxrel) CALL iso_verif_egalite_choix(xtent(iso_eau,il,im,jm), & qent(il,im,jm),'cv30_mixing 2112',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(xtelij(1,il,im,jm), & 'cv30_routine 2250') #endif enddo !do il = 1, nloc enddo !do jm = 1, klev enddo !do im = 1, klev #endif #endif #ifdef ISO #ifdef ISOTRAC ! seulement a la fin on taggue le condensat IF (option_cond.ge.1) THEN DO im = 1, nd DO jm = 1, nd DO il = 1, ncum ! colorier le condensat en un tag specifique DO ixt=niso+1,ntraciso IF (index_zone(ixt).EQ.izone_cond) THEN xtelij(ixt,il,im,jm)=xtelij(index_iso(ixt),il,im,jm) else !if (index_zone(ixt).EQ.izone_cond) THEN xtelij(ixt,il,im,jm)=0.0 endif !if (index_zone(ixt).EQ.izone_cond) THEN enddo !do ixt=1,ntraciso #ifdef ISOVERIF CALL iso_verif_egalite_choix(xtelij(iso_eau,il,im,jm), & elij(il,im,jm),'cv30_routines 2408',errmax,errmaxrel) CALL iso_verif_traceur(xtelij(1,il,im,jm), & 'condiso_liq_ice_vectiso_trac 358') #endif enddo !do il = 1, ncum enddo !do jm = 1, nd enddo !do im = 1, nd DO im = 1, nd DO il = 1, ncum ! colorier le condensat en un tag specifique DO ixt=niso+1,ntraciso IF (index_zone(ixt).EQ.izone_cond) THEN xtclw(ixt,il,im)=xtclw(index_iso(ixt),il,im) else !if (index_zone(ixt).EQ.izone_cond) THEN xtclw(ixt,il,im)=0.0 endif !if (index_zone(ixt).EQ.izone_cond) THEN enddo !do ixt=1,ntraciso #ifdef ISOVERIF CALL iso_verif_egalite_choix(xtclw(iso_eau,il,im), & clw(il,im),'cv30_routines 2427',errmax,errmaxrel) CALL iso_verif_traceur(xtclw(1,il,im), & 'condiso_liq_ice_vectiso_trac 358') IF (iso_verif_positif_nostop(xtclw(itZonIso( & izone_cond,iso_eau),i,k)-xtclw(iso_eau,i,k) & ,'cv30_routines 909').EQ.1) THEN WRITE(*,*) 'i,k=',i,k WRITE(*,*) 'xtclw=',xtclw(:,i,k) WRITE(*,*) 'niso,ntraciso,index_zone,izone_cond=', & niso,ntraciso,index_zone,izone_cond stop endif !if (iso_verif_positif_nostop(xtclw(itZonIso( #endif enddo !do il = 1, ncum enddo !do im = 1, nd ! WRITE(*,*) 'xtclw(:,1,2)=',xtclw(:,1,2) endif !if (option_tmin.EQ.1) THEN #endif #endif END SUBROUTINE cv30_mixing SUBROUTINE cv30_unsat(nloc, ncum, nd, na, ntra, icb, inb, t, rr, rs, gz, u, & v, tra, p, ph, th, tv, lv, cpn, ep, sigp, clw, m, ment, elij, delt, plcl, & mp, rp, up, vp, trap, wt, water, evap, b & ! RomP-jyg , wdtraina, wdtrainm & ! 26/08/10 RomP-jyg #ifdef ISO ,xt,xtclw,xtelij & ,xtp,xtwater,xtevap,xtwdtraina & #endif ) #ifdef ISO USE infotrac_phy, ONLY: ntraciso=>ntiso, niso USE isotopes_mod, ONLY: essai_convergence, iso_eau,iso_HDO,ridicule USE isotopes_routines_mod, ONLY: appel_stewart_vectall,appel_stewart_debug #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: errmax,errmaxrel, & iso_verif_egalite_choix, iso_verif_noNaN,iso_verif_aberrant, & iso_verif_egalite,iso_verif_egalite_choix_nostop,iso_verif_positif_nostop, & iso_verif_egalite_nostop,iso_verif_aberrant_nostop,deltaD,iso_verif_noNaN_nostop, & iso_verif_positif,iso_verif_egalite_vect2D #endif #ifdef ISOTRAC USE isotrac_mod, ONLY: option_cond,izone_cond USE infotrac_phy, ONLY: itZonIso #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: iso_verif_traceur_justmass, & & iso_verif_traceur USE isotrac_routines_mod, ONLY: iso_verif_traceur_pbidouille #endif #endif #endif USE cvflag_mod_h USE cvthermo_mod_h IMPLICIT NONE ! inputs: INTEGER ncum, nd, na, ntra, nloc INTEGER icb(nloc), inb(nloc) REAL delt, plcl(nloc) REAL t(nloc, nd), rr(nloc, nd), rs(nloc, nd) REAL u(nloc, nd), v(nloc, nd) REAL tra(nloc, nd, ntra) REAL p(nloc, nd), ph(nloc, nd+1) REAL th(nloc, na), gz(nloc, na) REAL lv(nloc, na), ep(nloc, na), sigp(nloc, na), clw(nloc, na) REAL cpn(nloc, na), tv(nloc, na) REAL m(nloc, na), ment(nloc, na, na), elij(nloc, na, na) #ifdef ISO REAL xt(ntraciso,nloc,nd), xtclw(ntraciso,nloc,na) REAL xtelij(ntraciso,nloc,na,na) ! real xtep(ntraciso,nloc,na) ! le 7 mai: on supprime xtep #endif ! outputs: REAL mp(nloc, na), rp(nloc, na), up(nloc, na), vp(nloc, na) REAL water(nloc, na), evap(nloc, na), wt(nloc, na) REAL trap(nloc, na, ntra) REAL b(nloc, na) ! 25/08/10 - RomP---- ajout des masses precipitantes ejectees ! lascendance adiabatique et des flux melanges Pa et Pm. ! Distinction des wdtrain ! Pa = wdtrainA Pm = wdtrainM REAL wdtraina(nloc, na), wdtrainm(nloc, na) #ifdef ISO REAL xtp(ntraciso,nloc,na) REAL xtwater(ntraciso,nloc,na), xtevap(ntraciso,nloc,na) REAL xtwdtraina(ntraciso,nloc,na) #endif ! local variables INTEGER i, j, k, il, num1 REAL tinv, delti REAL awat, afac, afac1, afac2, bfac REAL pr1, pr2, sigt, b6, c6, revap, tevap, delth REAL amfac, amp2, xf, tf, fac2, ur, sru, fac, d, af, bf REAL ampmax REAL lvcp(nloc, na) REAL wdtrain(nloc) LOGICAL lwork(nloc) #ifdef ISO INTEGER ixt REAL xtawat(ntraciso) REAL xtwdtrain(ntraciso,nloc) ! LOGICAL negation REAL rpprec(nloc,na) !#ifdef ISOVERIF ! integer iso_verif_aberrant_nostop !#ifdef ISOTRAC ! integer iso_verif_traceur_choix_nostop ! integer iso_verif_positif_nostop !#endif !#endif #endif ! ------------------------------------------------------ !#ifdef ISOVERIF ! WRITE(*,*) 'cv30_routines 2382: entree dans cv3_unsat' !#endif delti = 1./delt tinv = 1./3. mp(:, :) = 0. #ifdef ISO ! initialisation plus complete de water et rp water(:,:)=0.0 xtwater(:,:,:)=0.0 rp(:,:)=0.0 xtp(:,:,:)=0.0 #endif DO i = 1, nl DO il = 1, ncum mp(il, i) = 0.0 rp(il, i) = rr(il, i) up(il, i) = u(il, i) vp(il, i) = v(il, i) wt(il, i) = 0.001 water(il, i) = 0.0 evap(il, i) = 0.0 b(il, i) = 0.0 lvcp(il, i) = lv(il, i)/cpn(il, i) #ifdef ISO rpprec(il,i)=rp(il,i) DO ixt=1,ntraciso xtp(ixt,il,i)=xt(ixt,il,i) xtwater(ixt,il,i)=0.0 xtevap(ixt,il,i)=0.0 enddo !-- debug #ifdef ISOVERIF IF(iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xt(iso_eau,il,i),rr(il,i), & 'cv30_unsat 2245 ',errmax,errmaxrel) CALL iso_verif_egalite_choix(xtp(iso_eau,il,i),rp(il,i), & 'cv30_unsat 2247 ',errmax,errmaxrel) endif !IF(iso_eau.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur(xt(1,il,i),'cv30_routine 2410') CALL iso_verif_traceur(xtp(1,il,i),'cv30_routine 2411') #endif #endif #endif END DO END DO ! do k=1,ntra ! do i=1,nd ! do il=1,ncum ! trap(il,i,k)=tra(il,i,k) ! enddo ! enddo ! enddo ! RomP >>> DO i = 1, nd DO il = 1, ncum wdtraina(il, i) = 0.0 wdtrainm(il, i) = 0.0 END DO END DO ! RomP <<< ! *** check whether ep(inb)=0, if so, skip precipitating *** ! *** downdraft calculation *** DO il = 1, ncum lwork(il) = .TRUE. IF (ep(il,inb(il))<0.0001) lwork(il) = .FALSE. END DO CALL zilch(wdtrain, ncum) #ifdef ISO CALL zilch(xtwdtrain,ncum*ntraciso) #endif DO i = nl + 1, 1, -1 num1 = 0 DO il = 1, ncum IF (i<=inb(il) .AND. lwork(il)) num1 = num1 + 1 END DO IF (num1<=0) GO TO 400 ! *** integrate liquid water equation to find condensed water *** ! *** and condensed water flux *** ! *** begin downdraft loop *** ! *** calculate detrained precipitation *** DO il = 1, ncum IF (i<=inb(il) .AND. lwork(il)) THEN IF (cvflag_grav) THEN wdtrain(il) = grav*ep(il, i)*m(il, i)*clw(il, i) wdtraina(il, i) = wdtrain(il)/grav ! Pa 26/08/10 RomP #ifdef ISO DO ixt=1,ntraciso ! xtwdtrain(ixt,il)=grav*xtep(ixt,il,i)*m(il,i)*xtclw(ixt,il,i) xtwdtrain(ixt,il)=grav*ep(il,i)*m(il,i)*xtclw(ixt,il,i) enddo !--debug: #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xtwdtrain(iso_eau,il), & wdtrain(il),'cv30_routines 2313',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur(xtwdtrain(1,il),'cv30_routine 2480') #endif #endif !--end debug #endif ELSE wdtrain(il) = 10.0*ep(il, i)*m(il, i)*clw(il, i) wdtraina(il, i) = wdtrain(il)/10. ! Pa 26/08/10 RomP #ifdef ISO DO ixt=1,ntraciso ! xtwdtrain(ixt,il)=10.0*xtep(ixt,il,i)*m(il,i)*xtclw(ixt,il,i) xtwdtrain(ixt,il)=10.0*ep(il,i)*m(il,i)*xtclw(ixt,il,i) xtwdtraina(ixt,il, i) = xtwdtrain(ixt,il)/10. enddo #endif END IF END IF END DO IF (i>1) THEN DO j = 1, i - 1 DO il = 1, ncum IF (i<=inb(il) .AND. lwork(il)) THEN awat = elij(il, j, i) - (1.-ep(il,i))*clw(il, i) awat = amax1(awat, 0.0) #ifdef ISO ! precip mixed drafts computed from: xtawat/xtelij = awat/elij IF (elij(il,j,i).NE.0.0) THEN DO ixt=1,ntraciso xtawat(ixt)=xtelij(ixt,il,j,i)*(awat/elij(il,j,i)) xtawat(ixt)=amax1(xtawat(ixt),0.0) enddo !! xtawat(ixt)=amin1(xtawat(ixt),xtelij(ixt,il,j,i)) !security.. else DO ixt=1,ntraciso xtawat(ixt)=0.0 enddo !do ixt=1,niso endif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xtawat(iso_eau), & awat,'cv30_routines 2391',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur(xtawat(1),'cv30_routine 2522') #endif #endif #endif IF (cvflag_grav) THEN wdtrain(il) = wdtrain(il) + grav*awat*ment(il, j, i) #ifdef ISO DO ixt=1,ntraciso xtwdtrain(ixt,il)=xtwdtrain(ixt,il) & +grav*xtawat(ixt)*ment(il,j,i) enddo !do ixt=1,ntraciso #endif ELSE wdtrain(il) = wdtrain(il) + 10.0*awat*ment(il, j, i) #ifdef ISO DO ixt=1,ntraciso xtwdtrain(ixt,il)=xtwdtrain(ixt,il) & +10.0*xtawat(ixt)*ment(il,j,i) enddo !!do ixt=1,ntraciso #endif END IF !if (cvflag_grav) THEN #ifdef ISO !--debug: #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xtwdtrain(iso_eau,il), & wdtrain(il),'cv30_routines 2366',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur(xtwdtrain(1,il),'cv30_routine 2540') IF (option_cond.ge.1) THEN ! on verifie que tout le detrainement est tagge condensat IF (iso_verif_positif_nostop( & xtwdtrain(itZonIso(izone_cond,iso_eau),il) & -xtwdtrain(iso_eau,il), & 'cv30_routines 2795').EQ.1) THEN WRITE(*,*) 'xtwdtrain(:,il)=',xtwdtrain(:,il) WRITE(*,*) 'xtelij(:,il,j,i)=',xtelij(:,il,j,i) WRITE(*,*) 'xtclw(:,il,i)=',xtclw(:,il,i) stop endif !if (iso_verif_positif_nostop(Pxtisup(iso_eau,il)- endif !if (option_cond.ge.1) THEN #endif #endif #endif END IF !IF (i<=inb(il) .AND. lwork(il)) THEN END DO END DO DO il = 1, ncum IF (cvflag_grav) THEN wdtrainm(il, i) = wdtrain(il)/grav - wdtraina(il, i) ! Pm 26/08/10 RomP ELSE wdtrainm(il, i) = wdtrain(il)/10. - wdtraina(il, i) ! Pm 26/08/10 RomP END IF END DO END IF ! *** find rain water and evaporation using provisional *** ! *** estimates of rp(i)and rp(i-1) *** DO il = 1, ncum IF (i<=inb(il) .AND. lwork(il)) THEN wt(il, i) = 45.0 IF (i<inb(il)) THEN rp(il, i) = rp(il, i+1) + (cpd*(t(il,i+1)-t(il, & i))+gz(il,i+1)-gz(il,i))/lv(il, i) rp(il, i) = 0.5*(rp(il,i)+rr(il,i)) END IF rp(il, i) = amax1(rp(il,i), 0.0) rp(il, i) = amin1(rp(il,i), rs(il,i)) rp(il, inb(il)) = rr(il, inb(il)) IF (i==1) THEN afac = p(il, 1)*(rs(il,1)-rp(il,1))/(1.0E4+2000.0*p(il,1)*rs(il,1)) ELSE rp(il, i-1) = rp(il, i) + (cpd*(t(il,i)-t(il, & i-1))+gz(il,i)-gz(il,i-1))/lv(il, i) rp(il, i-1) = 0.5*(rp(il,i-1)+rr(il,i-1)) rp(il, i-1) = amin1(rp(il,i-1), rs(il,i-1)) rp(il, i-1) = amax1(rp(il,i-1), 0.0) afac1 = p(il, i)*(rs(il,i)-rp(il,i))/(1.0E4+2000.0*p(il,i)*rs(il,i) & ) afac2 = p(il, i-1)*(rs(il,i-1)-rp(il,i-1))/ & (1.0E4+2000.0*p(il,i-1)*rs(il,i-1)) afac = 0.5*(afac1+afac2) END IF IF (i==inb(il)) afac = 0.0 afac = amax1(afac, 0.0) bfac = 1./(sigd*wt(il,i)) ! jyg1 ! cc sigt=1.0 ! cc IF(i.ge.icb)sigt=sigp(i) ! prise en compte de la variation progressive de sigt dans ! les couches icb et icb-1: ! pour plcl<ph(i+1), pr1=0 & pr2=1 ! pour plcl>ph(i), pr1=1 & pr2=0 ! pour ph(i+1)<plcl<ph(i), pr1 est la proportion a cheval ! sur le nuage, et pr2 est la proportion sous la base du ! nuage. pr1 = (plcl(il)-ph(il,i+1))/(ph(il,i)-ph(il,i+1)) pr1 = max(0., min(1.,pr1)) pr2 = (ph(il,i)-plcl(il))/(ph(il,i)-ph(il,i+1)) pr2 = max(0., min(1.,pr2)) sigt = sigp(il, i)*pr1 + pr2 ! jyg2 b6 = bfac*50.*sigd*(ph(il,i)-ph(il,i+1))*sigt*afac c6 = water(il, i+1) + bfac*wdtrain(il) - 50.*sigd*bfac*(ph(il,i)-ph( & il,i+1))*evap(il, i+1) IF (c6>0.0) THEN revap = 0.5*(-b6+sqrt(b6*b6+4.*c6)) evap(il, i) = sigt*afac*revap water(il, i) = revap*revap ELSE evap(il, i) = -evap(il, i+1) + 0.02*(wdtrain(il)+sigd*wt(il,i)* & water(il,i+1))/(sigd*(ph(il,i)-ph(il,i+1))) END IF #ifdef ISO ! ajout cam: eviter les evaporations ou eaux negatives ! water(il,i)=max(0.0,water(il,i)) ! ceci est toujours verifie #ifdef ISOVERIF CALL iso_verif_positif(water(il,i),'cv30_unsat 2376') #endif ! evap(il,i)=max(0.0,evap(il,i)) ! evap<0 permet la conservation de ! l'eau ! fin ajout cam #endif ! *** calculate precipitating downdraft mass flux under *** ! *** hydrostatic approximation *** IF (i/=1) THEN tevap = amax1(0.0, evap(il,i)) delth = amax1(0.001, (th(il,i)-th(il,i-1))) IF (cvflag_grav) THEN mp(il, i) = 100.*ginv*lvcp(il, i)*sigd*tevap*(p(il,i-1)-p(il,i))/ & delth ELSE mp(il, i) = 10.*lvcp(il, i)*sigd*tevap*(p(il,i-1)-p(il,i))/delth END IF ! *** if hydrostatic assumption fails, *** ! *** solve cubic difference equation for downdraft theta *** ! *** and mass flux from two simultaneous differential eqns *** amfac = sigd*sigd*70.0*ph(il, i)*(p(il,i-1)-p(il,i))* & (th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i)) amp2 = abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i)) IF (amp2>(0.1*amfac)) THEN xf = 100.0*sigd*sigd*sigd*(ph(il,i)-ph(il,i+1)) tf = b(il, i) - 5.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i)* & sigd*th(il,i)) af = xf*tf + mp(il, i+1)*mp(il, i+1)*tinv bf = 2.*(tinv*mp(il,i+1))**3 + tinv*mp(il, i+1)*xf*tf + & 50.*(p(il,i-1)-p(il,i))*xf*tevap fac2 = 1.0 IF (bf<0.0) fac2 = -1.0 bf = abs(bf) ur = 0.25*bf*bf - af*af*af*tinv*tinv*tinv IF (ur>=0.0) THEN sru = sqrt(ur) fac = 1.0 IF ((0.5*bf-sru)<0.0) fac = -1.0 mp(il, i) = mp(il, i+1)*tinv + (0.5*bf+sru)**tinv + & fac*(abs(0.5*bf-sru))**tinv ELSE d = atan(2.*sqrt(-ur)/(bf+1.0E-28)) IF (fac2<0.0) d = 3.14159 - d mp(il, i) = mp(il, i+1)*tinv + 2.*sqrt(af*tinv)*cos(d*tinv) END IF mp(il, i) = amax1(0.0, mp(il,i)) IF (cvflag_grav) THEN ! jyg : il y a vraisemblablement une erreur dans la ligne 2 ! suivante: ! il faut diviser par (mp(il,i)*sigd*grav) et non par ! (mp(il,i)+sigd*0.1). ! Et il faut bien revoir les facteurs 100. b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap/(mp(il, & i)+sigd*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i & )*sigd*th(il,i)) ELSE b(il, i-1) = b(il, i) + 100.0*(p(il,i-1)-p(il,i))*tevap/(mp(il, & i)+sigd*0.1) - 10.0*(th(il,i)-th(il,i-1))*t(il, i)/(lvcp(il,i & )*sigd*th(il,i)) END IF b(il, i-1) = amax1(b(il,i-1), 0.0) END IF ! *** limit magnitude of mp(i) to meet cfl condition ! *** ampmax = 2.0*(ph(il,i)-ph(il,i+1))*delti amp2 = 2.0*(ph(il,i-1)-ph(il,i))*delti ampmax = amin1(ampmax, amp2) mp(il, i) = amin1(mp(il,i), ampmax) ! *** force mp to decrease linearly to zero ! *** ! *** between cloud base and the surface ! *** IF (p(il,i)>p(il,icb(il))) THEN mp(il, i) = mp(il, icb(il))*(p(il,1)-p(il,i))/ & (p(il,1)-p(il,icb(il))) END IF END IF ! i.EQ.1 ! *** find mixing ratio of precipitating downdraft *** IF (i/=inb(il)) THEN rp(il, i) = rr(il, i) IF (mp(il,i)>mp(il,i+1)) THEN IF (cvflag_grav) THEN rp(il, i) = rp(il, i+1)*mp(il, i+1) + & rr(il, i)*(mp(il,i)-mp(il,i+1)) + 100.*ginv*0.5*sigd*(ph(il,i & )-ph(il,i+1))*(evap(il,i+1)+evap(il,i)) ELSE rp(il, i) = rp(il, i+1)*mp(il, i+1) + & rr(il, i)*(mp(il,i)-mp(il,i+1)) + 5.*sigd*(ph(il,i)-ph(il,i+1 & ))*(evap(il,i+1)+evap(il,i)) END IF rp(il, i) = rp(il, i)/mp(il, i) up(il, i) = up(il, i+1)*mp(il, i+1) + u(il, i)*(mp(il,i)-mp(il,i+ & 1)) up(il, i) = up(il, i)/mp(il, i) vp(il, i) = vp(il, i+1)*mp(il, i+1) + v(il, i)*(mp(il,i)-mp(il,i+ & 1)) vp(il, i) = vp(il, i)/mp(il, i) ! do j=1,ntra ! trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1) ! testmaf : +trap(il,i,j)*(mp(il,i)-mp(il,i+1)) ! : +tra(il,i,j)*(mp(il,i)-mp(il,i+1)) ! trap(il,i,j)=trap(il,i,j)/mp(il,i) ! END DO ELSE IF (mp(il,i+1)>1.0E-16) THEN IF (cvflag_grav) THEN rp(il, i) = rp(il, i+1) + 100.*ginv*0.5*sigd*(ph(il,i)-ph(il, & i+1))*(evap(il,i+1)+evap(il,i))/mp(il, i+1) ELSE rp(il, i) = rp(il, i+1) + 5.*sigd*(ph(il,i)-ph(il,i+1))*(evap & (il,i+1)+evap(il,i))/mp(il, i+1) END IF up(il, i) = up(il, i+1) vp(il, i) = vp(il, i+1) ! do j=1,ntra ! trap(il,i,j)=trap(il,i+1,j) ! END DO END IF END IF #ifdef ISO rpprec(il,i)=max(rp(il,i),0.0) #endif rp(il, i) = amin1(rp(il,i), rs(il,i)) rp(il, i) = amax1(rp(il,i), 0.0) END IF END IF END DO #ifdef ISO #ifdef ISOVERIF ! verif des inputs a appel stewart ! WRITE(*,*) 'cv30_routines 2842 tmp: appel de appel_stewart' DO il=1,ncum IF (i.le.inb(il) .AND. lwork(il)) THEN IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xt(iso_eau,il,i), & rr(il,i),'appel_stewart 262, cas 1.1',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN !#ifdef ISOTRAC ! if (option_tmin.ge.1) THEN ! CALL iso_verif_positif(xtwater( ! : itZonIso(izone_cond,iso_eau),il,i+1) ! : -xtwater(iso_eau,il,i+1), ! : 'cv30_routines 3083') ! endif !if (option_tmin.ge.1) THEN !#endif endif enddo #endif IF (1.EQ.0) THEN ! appel de appel_stewart_vectorise CALL appel_stewart_vectall(lwork,ncum, & ph,t,evap,xtwdtrain, & wdtrain, & water,rr,xt,rs,rpprec,mp,wt, & ! inputs physiques xtwater,xtp, & ! outputs indispensables xtevap, & ! diagnostiques sigd, & ! inputs tunables i,inb, & ! altitude: car cas particulier en INB na,nd,nloc,cvflag_grav,ginv,1e-16) else !if (1.EQ.0) THEN ! truc simple sans fractionnement ! juste pour debuggage CALL appel_stewart_debug(lwork,nloc,inb,na,i, & evap,water,rpprec,rr,wdtrain, & xtevap,xtwater,xtp,xt,xtwdtrain) endif ! if (1.EQ.0) THEN #ifdef ISOVERIF ! WRITE(*,*) 'cv30_routines 2864 tmp: sortie de appel_stewart' ! verif des outputs de appel stewart DO il=1,ncum IF (i.le.inb(il) .AND. lwork(il)) THEN DO ixt=1,ntraciso CALL iso_verif_noNAN(xtp(ixt,il,i),'cv30_unsat 3382') CALL iso_verif_noNAN(xtwater(ixt,il,i),'cv30_unsat 3381') CALL iso_verif_noNAN(xtevap(ixt,il,i),'cv30_unsat 2661') enddo IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xtp(iso_eau,il,i), & rpprec(il,i),'cv30_unsat 2736',errmax,errmaxrel) CALL iso_verif_egalite_choix(xtwater(iso_eau,il,i), & water(il,i),'cv30_unsat 2747',errmax,errmaxrel) ! WRITE(*,*) 'xtwater(4,il,i)=',xtwater(4,il,i) ! WRITE(*,*) 'water(il,i)=',water(il,i) CALL iso_verif_egalite_choix(xtevap(iso_eau,il,i), & evap(il,i),'cv30_unsat 2751',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN IF ((iso_HDO.gt.0).AND. & (rp(il,i).gt.ridicule)) THEN CALL iso_verif_aberrant(xtp(iso_HDO,il,i)/rpprec(il,i), & 'cv3unsat 2756') endif !if ((iso_HDO.gt.0).AND. #ifdef ISOTRAC ! if (il.EQ.602) THEN ! WRITE(*,*) 'cv30_routine tmp: il,i=',il,i ! WRITE(*,*) 'xtp(iso_eau:ntraciso:3,il,i)=', ! : xtp(iso_eau:ntraciso:3,il,i) ! endif CALL iso_verif_traceur(xtp(1,il,i),'cv30_routine 2852') CALL iso_verif_traceur(xtwater(1,il,1), & 'cv30_routine 2853 unsat apres appel') CALL iso_verif_traceur_pbidouille(xtwater(1,il,i), & 'cv30_routine 2853b') CALL iso_verif_traceur_justmass(xtevap(1,il,i), & 'cv30_routine 2854') ! if (option_tmin.ge.1) THEN ! CALL iso_verif_positif(xtwater( ! : itZonIso(izone_cond,iso_eau),il,i) ! : -xtwater(iso_eau,il,i), ! : 'cv30_routines 3143') ! endif !if (option_tmin.ge.1) THEN #endif endif !if (i.le.inb(il) .AND. lwork(il)) THEN enddo !do il=1,ncum #endif ! equivalent isotopique de rp(il,i)=amin1(rp(il,i),rs(il,i)) DO il=1,ncum IF (i.lt.inb(il) .AND. lwork(il)) THEN IF (rpprec(il,i).gt.rs(il,i)) THEN IF (rs(il,i).le.0) THEN WRITE(*,*) 'cv3unsat 2640' stop endif DO ixt=1,ntraciso xtp(ixt,il,i)=xtp(ixt,il,i)/rpprec(il,i)*rs(il,i) xtp(ixt,il,i)=max(0.0,xtp(ixt,il,i)) enddo !do ixt=1,niso #ifdef ISOVERIF DO ixt=1,ntraciso CALL iso_verif_noNaN(xtp(ixt,il,i),'cv3unsat 2641') enddo !do ixt=1,niso IF (iso_eau.gt.0) THEN ! WRITE(*,*) 'xtp(iso_eau,il,i)=',xtp(iso_eau,il,i) CALL iso_verif_egalite_choix(xtp(iso_eau,il,i),rp(il,i), & 'cv3unsat 2653',errmax,errmaxrel) CALL iso_verif_egalite_choix(xtp(iso_eau,il,i), & rs(il,i),'cv3unsat 2654',errmax,errmaxrel) endif IF ((iso_HDO.gt.0).AND. & (rp(il,i).gt.ridicule)) THEN IF (iso_verif_aberrant_nostop(xtp(iso_HDO,il,i)/rp(il,i), & 'cv3unsat 2658').EQ.1) THEN WRITE(*,*) 'rpprec(il,i),rs(il,i),rp(il,i)=', & rpprec(il,i),rs(il,i),rp(il,i) stop endif endif #ifdef ISOTRAC CALL iso_verif_traceur(xtp(1,il,i),'cv30_routine 2893') #endif #endif rpprec(il,i)=rs(il,i) endif !if (rp(il,i).gt.rs(il,i)) THEN endif !if (i.lt.INB et lwork) enddo ! il=1,ncum #endif 400 END DO ! fin de la boucle en i (altitude) #ifdef ISO WRITE(*,*) 'nl=',nl,'nd=',nd,'; ncum=',ncum #ifdef ISOVERIF DO i=1,nl !nl DO il=1,ncum IF (iso_eau.gt.0) THEN ! WRITE(*,*) 'cv30_routines 2767:i,il,lwork(il),inb(il)=', ! : i,il,lwork(il),inb(il) ! WRITE(*,*) 'rp(il,i),xtp(iso_eau,il,i)=', ! : rp(il,i),xtp(iso_eau,il,i) CALL iso_verif_egalite_choix(xt(iso_eau,il,i), & rr(il,i),'cv30_unsat 2668',errmax,errmaxrel) CALL iso_verif_egalite_choix(xtp(iso_eau,il,i), & rp(il,i),'cv30_unsat 2670',errmax,errmaxrel) CALL iso_verif_egalite_choix(xtwater(iso_eau,il,i), & water(il,i),'cv30_unsat 2672',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN !#ifdef ISOTRAC ! if (iso_verif_traceur_choix_nostop(xtwater(1,il,i), ! : 'cv30_routine 2982 unsat',errmax, ! : errmaxrel,ridicule_trac,deltalimtrac).EQ.1) THEN ! WRITE(*,*) 'il,i,inb(il),lwork(il)=', ! : il,i,inb(il),lwork(il) ! WRITE(*,*) 'xtwater(:,il,i)=',xtwater(:,il,i) ! stop ! endif !#endif enddo !do il=1,nloc!ncum enddo !do i=1,nl!nl il=5 i=39 WRITE(*,*) 'cv30_unsat 2780: il,water(il,i),xtwater(iso_eau,il,i)=' & ,il,water(il,i),xtwater(iso_eau,il,i) #endif #endif END SUBROUTINE cv30_unsat SUBROUTINE cv30_yield(nloc, ncum, nd, na, ntra, icb, inb, delt, t, rr, u, v, & tra, gz, p, ph, h, hp, lv, cpn, th, ep, clw, m, tp, mp, rp, up, vp, trap, & wt, water, evap, b, ment, qent, uent, vent, nent, elij, traent, sig, tv, & tvp, iflag, precip, vprecip, ft, fr, fu, fv, ftra, upwd, dnwd, dnwd0, ma, & mike, tls, tps, qcondc, wd & #ifdef ISO ,xt,xtclw,xtp,xtwater,xtevap & ,xtent,xtelij,xtprecip,fxt,xtVprecip & #ifdef DIAGISO ,fq_detrainement,fq_ddft,fq_fluxmasse,fq_evapprecip & ,fxt_detrainement,fxt_ddft,fxt_fluxmasse,fxt_evapprecip & ,f_detrainement,q_detrainement,xt_detrainement & #endif #endif ) USE conema3_mod_h #ifdef ISO USE infotrac_phy, ONLY: ntraciso=>ntiso, niso, nzone, itZonIso USE isotopes_mod, ONLY: essai_convergence,ridicule,iso_eau,iso_HDO,iso_O18 #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: errmax,errmaxrel, & iso_verif_egalite_choix, iso_verif_noNaN,iso_verif_aberrant, & iso_verif_egalite,iso_verif_egalite_choix_nostop,iso_verif_positif_nostop, & iso_verif_egalite_nostop,iso_verif_aberrant_nostop,deltaD,iso_verif_noNaN_nostop, & iso_verif_positif,iso_verif_egalite_vect2D, & iso_verif_aberrant_enc_nostop,iso_verif_aberrant_encadre,iso_verif_o18_aberrant, & iso_verif_O18_aberrant_nostop,deltaO #endif #ifdef ISOTRAC USE isotrac_mod, ONLY: option_traceurs, & izone_revap,izone_poubelle,izone_ddft #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: iso_verif_traceur_choix_nostop,deltalimtrac, & & iso_verif_tracpos_choix_nostop,iso_verif_traceur,iso_verif_traceur_justmass USE isotrac_mod, ONLY: ridicule_trac #endif #endif #endif USE cvflag_mod_h USE cvthermo_mod_h IMPLICIT NONE ! inputs: INTEGER ncum, nd, na, ntra, nloc INTEGER icb(nloc), inb(nloc) REAL delt REAL t(nloc, nd), rr(nloc, nd), u(nloc, nd), v(nloc, nd) REAL tra(nloc, nd, ntra), sig(nloc, nd) REAL gz(nloc, na), ph(nloc, nd+1), h(nloc, na), hp(nloc, na) REAL th(nloc, na), p(nloc, nd), tp(nloc, na) REAL lv(nloc, na), cpn(nloc, na), ep(nloc, na), clw(nloc, na) REAL m(nloc, na), mp(nloc, na), rp(nloc, na), up(nloc, na) REAL vp(nloc, na), wt(nloc, nd), trap(nloc, nd, ntra) REAL water(nloc, na), evap(nloc, na), b(nloc, na) REAL ment(nloc, na, na), qent(nloc, na, na), uent(nloc, na, na) ! ym real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na) REAL vent(nloc, na, na), elij(nloc, na, na) INTEGER nent(nloc, na) REAL traent(nloc, na, na, ntra) REAL tv(nloc, nd), tvp(nloc, nd) #ifdef ISO REAL xt(ntraciso,nloc,nd) ! real xtep(ntraciso,nloc,na) ! le 7 mai: on supprime xtep REAL xtclw(ntraciso,nloc,na), xtp(ntraciso,nloc,na) REAL xtwater(ntraciso,nloc,na), xtevap(ntraciso,nloc,na) REAL xtent(ntraciso,nloc,na,na), xtelij(ntraciso,nloc,na,na) #ifdef ISOVERIF CHARACTER (LEN=20) :: modname='cv30_compress' CHARACTER (LEN=80) :: abort_message #endif #endif ! input/output: INTEGER iflag(nloc) ! outputs: REAL precip(nloc) REAL vprecip(nloc, nd+1) REAL ft(nloc, nd), fr(nloc, nd), fu(nloc, nd), fv(nloc, nd) REAL ftra(nloc, nd, ntra) REAL upwd(nloc, nd), dnwd(nloc, nd), ma(nloc, nd) REAL dnwd0(nloc, nd), mike(nloc, nd) REAL tls(nloc, nd), tps(nloc, nd) REAL qcondc(nloc, nd) ! cld REAL wd(nloc) ! gust #ifdef ISO REAL xtprecip(ntraciso,nloc), fxt(ntraciso,nloc,nd) REAL xtVprecip(ntraciso,nloc,nd+1) #endif ! local variables: INTEGER i, k, il, n, j, num1 REAL rat, awat, delti REAL ax, bx, cx, dx, ex REAL cpinv, rdcp, dpinv REAL lvcp(nloc, na), mke(nloc, na) REAL am(nloc), work(nloc), ad(nloc), amp1(nloc) ! !! real up1(nloc), dn1(nloc) REAL up1(nloc, nd, nd), dn1(nloc, nd, nd) REAL asum(nloc), bsum(nloc), csum(nloc), dsum(nloc) REAL qcond(nloc, nd), nqcond(nloc, nd), wa(nloc, nd) ! cld REAL siga(nloc, nd), sax(nloc, nd), mac(nloc, nd) ! cld #ifdef ISO INTEGER ixt REAL xtbx(ntraciso), xtawat(ntraciso) ! cam debug ! pour l'homogeneisation sous le nuage: REAL frsum(nloc), bxtsum(ntraciso,nloc), fxtsum(ntraciso,nloc) ! correction dans calcul tendance liee a Am: REAL dq_tmp,k_tmp,dx_tmp,R_tmp,dqreste_tmp,dxreste_tmp,kad_tmp LOGICAL correction_excess_aberrant parameter (correction_excess_aberrant=.FALSE.) ! correction qui permettait d'eviter deltas et dexcess aberrants. Mais ! pb: ne conserve pas la masse d'isotopes! #ifdef DIAGISO ! diagnostiques juste: tendance des differents processus REAL fxt_detrainement(ntraciso,nloc,nd) REAL fxt_fluxmasse(ntraciso,nloc,nd) REAL fxt_evapprecip(ntraciso,nloc,nd) REAL fxt_ddft(ntraciso,nloc,nd) REAL fq_detrainement(nloc,nd) REAL q_detrainement(nloc,nd) REAL xt_detrainement(ntraciso,nloc,nd) REAL f_detrainement(nloc,nd) REAL fq_fluxmasse(nloc,nd) REAL fq_evapprecip(nloc,nd) REAL fq_ddft(nloc,nd) #endif !#ifdef ISOVERIF ! integer iso_verif_aberrant_nostop ! real deltaD !#endif #ifdef ISOTRAC ! integer iso_verif_traceur_choix_nostop ! integer iso_verif_tracpos_choix_nostop REAL xtnew(ntraciso) ! real conversion(niso) REAL fxtYe(niso) REAL fxtqe(niso) REAL fxtXe(niso) REAL fxt_revap(niso) REAL Xe(niso) INTEGER ixt_revap,izone INTEGER ixt_poubelle, ixt_ddft,iiso #endif #endif ! ------------------------------------------------------------- ! initialization: delti = 1.0/delt DO il = 1, ncum precip(il) = 0.0 wd(il) = 0.0 ! gust vprecip(il, nd+1) = 0. #ifdef ISO ! cam debug ! WRITE(*,*) 'cv30_routines 3082: entree dans cv3_yield' ! en cam debug DO ixt = 1, ntraciso xtprecip(ixt,il)=0.0 xtVprecip(ixt,il,nd+1)=0.0 enddo #endif END DO DO i = 1, nd DO il = 1, ncum vprecip(il, i) = 0.0 ft(il, i) = 0.0 fr(il, i) = 0.0 fu(il, i) = 0.0 fv(il, i) = 0.0 qcondc(il, i) = 0.0 ! cld qcond(il, i) = 0.0 ! cld nqcond(il, i) = 0.0 ! cld #ifdef ISO DO ixt = 1, ntraciso fxt(ixt,il,i)=0.0 xtVprecip(ixt,il,i)=0.0 enddo #ifdef DIAGISO fq_fluxmasse(il,i)=0.0 fq_detrainement(il,i)=0.0 f_detrainement(il,i)=0.0 q_detrainement(il,i)=0.0 fq_evapprecip(il,i)=0.0 fq_ddft(il,i)=0.0 DO ixt = 1, niso fxt_fluxmasse(ixt,il,i)=0.0 fxt_detrainement(ixt,il,i)=0.0 xt_detrainement(ixt,il,i)=0.0 fxt_evapprecip(ixt,il,i)=0.0 fxt_ddft(ixt,il,i)=0.0 enddo #endif #endif END DO END DO ! do j=1,ntra ! do i=1,nd ! do il=1,ncum ! ftra(il,i,j)=0.0 ! enddo ! enddo ! enddo DO i = 1, nl DO il = 1, ncum lvcp(il, i) = lv(il, i)/cpn(il, i) END DO END DO ! *** calculate surface precipitation in mm/day *** DO il = 1, ncum IF (ep(il,inb(il))>=0.0001) THEN IF (cvflag_grav) THEN precip(il) = wt(il, 1)*sigd*water(il, 1)*86400.*1000./(rowl*grav) #ifdef ISO DO ixt = 1, ntraciso xtprecip(ixt,il)=wt(il,1)*sigd*xtwater(ixt,il,1) & *86400.*1000./(rowl*grav) ! en mm/jour enddo ! cam verif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN ! WRITE(*,*) 'cv30_yield 2952: '// ! : 'il,water(il,1),xtwater(iso_eau,il,1)=' ! : ,il,water(il,1),xtwater(iso_eau,il,1) CALL iso_verif_egalite_choix(xtwater(iso_eau,il,1), & water(il,1),'cv30_routines 2959', & errmax,errmaxrel) !Rq: wt(il,1)*sigd*86400.*1000./(rowl*grav)=3964.6565 ! -> on auatorise 3e3 fois plus d'erreur dans precip CALL iso_verif_egalite_choix(xtprecip(iso_eau,il), & precip(il),'cv30_routines 3138', & errmax*4e3,errmaxrel) endif !if (iso_eau.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur(xtwater(1,il,1), & 'cv30_routine 3146') IF (iso_verif_traceur_choix_nostop(xtprecip(1,il), & 'cv30_routine 3147',errmax*1e2, & errmaxrel,ridicule_trac,deltalimtrac).EQ.1) THEN WRITE(*,*) 'il,inb(il)=',il,inb(il) WRITE(*,*) 'xtwater(:,il,1)=',xtwater(:,il,1) WRITE(*,*) 'xtprecip(:,il)=',xtprecip(:,il) WRITE(*,*) 'fac=',wt(il,1)*sigd*86400.*1000./(rowl*grav) stop endif #endif #endif ! end cam verif #endif ELSE precip(il) = wt(il, 1)*sigd*water(il, 1)*8640. #ifdef ISO DO ixt = 1, ntraciso xtprecip(ixt,il)=wt(il,1)*sigd*xtwater(ixt,il,1)*8640. enddo ! cam verif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xtprecip(iso_eau,il), & precip(il),'cv30_routines 3139', & errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur(xtprecip(1,il),'cv30_routine 3166') #endif #endif ! end cam verif #endif END IF !IF (cvflag_grav) THEN END IF !IF (cvflag_grav) THEN END DO ! *** CALCULATE VERTICAL PROFILE OF PRECIPITATIONs IN kg/m2/s === ! MAF rajout pour lessivage DO k = 1, nl DO il = 1, ncum IF (k<=inb(il)) THEN IF (cvflag_grav) THEN vprecip(il, k) = wt(il, k)*sigd*water(il, k)/grav #ifdef ISO DO ixt=1,ntraciso xtVPrecip(ixt,il,k) = wt(il,k)*sigd & *xtwater(ixt,il,k)/grav enddo #endif ELSE vprecip(il, k) = wt(il, k)*sigd*water(il, k)/10. #ifdef ISO DO ixt=1,ntraciso xtVPrecip(ixt,il,k) = wt(il,k)*sigd & *xtwater(ixt,il,k)/10.0 enddo #endif END IF END IF END DO END DO ! *** Calculate downdraft velocity scale *** ! *** NE PAS UTILISER POUR L'INSTANT *** ! do il=1,ncum ! wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il)) ! : /(sigd*p(il,icb(il))) ! enddo ! *** calculate tendencies of lowest level potential temperature *** ! *** and mixing ratio *** DO il = 1, ncum work(il) = 1.0/(ph(il,1)-ph(il,2)) am(il) = 0.0 END DO DO k = 2, nl DO il = 1, ncum IF (k<=inb(il)) THEN am(il) = am(il) + m(il, k) END IF END DO END DO DO il = 1, ncum ! convect3 if((0.1*dpinv*am).ge.delti)iflag(il)=4 IF (cvflag_grav) THEN IF ((0.01*grav*work(il)*am(il))>=delti) iflag(il) = 1 !consist vect ft(il, 1) = 0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1)+(gz(il,2)-gz(il, & 1))/cpn(il,1)) ELSE IF ((0.1*work(il)*am(il))>=delti) iflag(il) = 1 !consistency vect ft(il, 1) = 0.1*work(il)*am(il)*(t(il,2)-t(il,1)+(gz(il,2)-gz(il, & 1))/cpn(il,1)) END IF ft(il, 1) = ft(il, 1) - 0.5*lvcp(il, 1)*sigd*(evap(il,1)+evap(il,2)) IF (cvflag_grav) THEN ft(il, 1) = ft(il, 1) - 0.009*grav*sigd*mp(il, 2)*t(il, 1)*b(il, 1)* & work(il) ELSE ft(il, 1) = ft(il, 1) - 0.09*sigd*mp(il, 2)*t(il, 1)*b(il, 1)*work(il) END IF ft(il, 1) = ft(il, 1) + 0.01*sigd*wt(il, 1)*(cl-cpd)*water(il, 2)*(t(il,2 & )-t(il,1))*work(il)/cpn(il, 1) IF (cvflag_grav) THEN ! jyg1 Correction pour mieux conserver l'eau (conformite avec ! CONVECT4.3) ! (sb: pour l'instant, on ne fait que le chgt concernant grav, pas ! evap) fr(il, 1) = 0.01*grav*mp(il, 2)*(rp(il,2)-rr(il,1))*work(il) + & sigd*0.5*(evap(il,1)+evap(il,2)) ! +tard : +sigd*evap(il,1) fr(il, 1) = fr(il, 1) + 0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il) #ifdef ISO ! juste Mp et evap pour l'instant, voir plus bas pour am DO ixt = 1, ntraciso fxt(ixt,il,1)= & 0.01*grav*mp(il,2)*(xtp(ixt,il,2)-xt(ixt,il,1))*work(il) & +sigd*0.5*(xtevap(ixt,il,1)+xtevap(ixt,il,2)) !c+tard : +sigd*xtevap(ixt,il,1) enddo !do ixt = 1, ntraciso ! pour water tagging option 6: pas besoin ici de faire de conversion. #ifdef DIAGISO fq_ddft(il,1)=fq_ddft(il,1) & +0.01*grav*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) fq_evapprecip(il,1)=fq_evapprecip(il,1) & +sigd*0.5*(evap(il,1)+evap(il,2)) fq_fluxmasse(il,1)=fq_fluxmasse(il,1) & +0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il) DO ixt = 1, ntraciso ! fxt_fluxmasse(ixt,il,1)=fxt_fluxmasse(ixt,il,1) & ! & +0.01*grav*am(il)*(xt(ixt,il,2)-xt(ixt,il,1))*work(il) ! deplace ! plus haut car il existe differents cas fxt_ddft(ixt,il,1)=fxt_ddft(ixt,il,1) & +0.01*grav*mp(il,2)*(xtp(ixt,il,2)-xt(ixt,il,1))*work(il) fxt_evapprecip(ixt,il,1)=fxt_evapprecip(ixt,il,1) & +sigd*0.5*(xtevap(ixt,il,1)+xtevap(ixt,il,2)) enddo #endif ! pour l'ajout de la tendance liee au flux de masse Am, il faut etre ! prudent. ! On a dq1=k*(q2-q1) avec k=dt*0.01*grav*am(il)*work(il) ! Pour les isotopes, la formule utilisee depuis 2006 et qui avait toujours marche est: ! dx1=k*(x2-x1) ! Mais on plante dans un cas pathologique en decembre 2017 lors du test ! d'un cas d'Anne Cozic: les isotopes deviennent negatifs. ! C'est un cas pas physique: on perd 99% de la masse de vapeur d'eau! ! q2=1.01e-3 et q1=1.25e-3 kg/kg ! et dq=-1.24e-3: comment est-ce possible qu'un flux venant d'un air a ! q2= 1.01e-3 asseche q1 jusqu'a 0.01e-3kg/kg! ! Pour les isotopes, ca donne des x1+dx negatifs. ! Ce n'est pas physique mais il faut quand meme s'adapter. ! Pour cela, on considere que d'abord on fait rentrer le flux de masse ! descendant, et ensuite seulement on fait sortir le flux de masse ! sortant. ! Ainsi, le flux de masse sortant ne modifie pas la composition ! isotopique de la vapeur d'eau q1. ! A la fin, on a R=(x1+dx)/(q1+dq)=(x1+k*x2)/(q1+k*q2) ! On verifie que quand k est petit, on tend vers la formulation ! habituelle. ! Comme on est habitues a la formulation habituelle, qu'elle a fait ses ! preuves, on la garde sauf dans le cas ou dq/q<-0.9 ou on utilise la ! nouvelle formulation. ! rappel: dq_tmp=0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)*delt ! Meme avec cette nouvelle foirmulation, on a encore des isotopes ! negatifs, cette fois a cause des ddfts ! On considere donc les tendances et serie et non en parallele quand on ! calcule R_tmp. dq_tmp=0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)*delt ! utile ci-dessous IF ((dq_tmp/rr(il,1).lt.-0.9).AND.correction_excess_aberrant) THEN ! nouvelle formulation ou on fait d'abord entrer k*q2 et ensuite ! seulement on fait sortir k*q1 sans changement de composition ! isotopique k_tmp=0.01*grav*am(il)*work(il)*delt dqreste_tmp= 0.01*grav*mp(il, 2)*(rp(il,2)-rr(il,1))*work(il)*delt + & sigd*0.5*(evap(il,1)+evap(il,2))*delt DO ixt = 1, ntraciso dxreste_tmp= 0.01*grav*mp(il,2)*(xtp(ixt,il,2)-xt(ixt,il,1))*work(il)*delt & +sigd*0.5*(xtevap(ixt,il,1)+xtevap(ixt,il,2))*delt R_tmp=(xt(ixt,il,1)+dxreste_tmp+k_tmp*xt(ixt,il,2))/(rr(il,1)+dqreste_tmp+k_tmp*rr(il,2)) dx_tmp=R_tmp*(rr(il,1)+dqreste_tmp+dq_tmp)-(xt(ixt,il,1)+dxreste_tmp) fxt(ixt,il,1)=fxt(ixt,il,1) & + dx_tmp/delt #ifdef ISOVERIF IF (ixt.EQ.iso_HDO) THEN WRITE(*,*) 'cv30_routines 3888: il=',il WRITE(*,*) 'dq_tmp,rr(il,1)=',dq_tmp,rr(il,1) WRITE(*,*) 'R_tmp,dx_tmp,delt=',R_tmp,dx_tmp,delt WRITE(*,*) 'xt(ixt,il,1:2)=',xt(ixt,il,1:2) WRITE(*,*) 'rr(il,1:2)=',rr(il,1:2) WRITE(*,*) 'fxt=',dx_tmp/delt WRITE(*,*) 'rr(il,1)+dq_tmp=',rr(il,1)+dq_tmp WRITE(*,*) 'xt(ixt,il,1)+dx_tmp=',xt(ixt,il,1)+dx_tmp WRITE(*,*) 'xt(ixt,il,1)+fxt(ixt,il,1)*delt=', & xt(ixt,il,1)+fxt(ixt,il,1)*delt WRITE(*,*) 'dqreste_tmp,dxreste_tmp=',dqreste_tmp,dxreste_tmp WRITE(*,*) 'formule classique: fxt_Am=',0.01*grav*am(il)*(xt(ixt,il,2)-xt(ixt,il,1))*work(il) WRITE(*,*) 'donnerait dxt=',0.01*grav*am(il)*(xt(ixt,il,2)-xt(ixt,il,1))*work(il)*delt endif !if (ixt.EQ.iso_HDO) THEN #endif #ifdef DIAGISO IF (ixt.le.niso) THEN fxt_fluxmasse(ixt,il,1)=fxt_fluxmasse(ixt,il,1) & + dx_tmp/delt endif #endif enddo ! do ixt = 1, ntraciso else !if (dq_tmp/rr(il,1).lt.-0.9) THEN ! formulation habituelle qui avait toujours marche de 2006 a ! decembre 2017. DO ixt = 1, ntraciso fxt(ixt,il,1)=fxt(ixt,il,1) & +0.01*grav*am(il)*(xt(ixt,il,2)-xt(ixt,il,1))*work(il) #ifdef DIAGISO IF (ixt.le.niso) THEN fxt_fluxmasse(ixt,il,1)=fxt_fluxmasse(ixt,il,1) & +0.01*grav*am(il)*(xt(ixt,il,2)-xt(ixt,il,1))*work(il) endif #endif enddo !do ixt = 1, ntraciso endif !if (dq_tmp/rr(il,1).lt.-0.9) THEN ! cam verif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,1), & fr(il,1),'cv30_routines 3251', & errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN !WRITE(*,*) 'il,am(il)=',il,am(il) IF ((iso_HDO.gt.0).AND. & (rr(il,1)+delt*fr(il,1).gt.ridicule)) THEN IF (iso_verif_aberrant_enc_nostop((xt(iso_HDO,il,1) & +delt*fxt(iso_HDO,il,1))/(rr(il,1)+delt*fr(il,1)), & 'cv30_yield 3125, ddft en 1').EQ.1) THEN WRITE(*,*) 'il,rr(il,1),delt=',il,rr(il,1),delt WRITE(*,*) 'deltaDxt=',deltaD(xt(iso_HDO,il,1)/rr(il,1)) WRITE(*,*) 'delt*fr(il,1),fr(il,1)=',delt*fr(il,1),fr(il,1) WRITE(*,*) 'fxt=',fxt(iso_HDO,il,1) WRITE(*,*) 'fq_ddft(il,1)=',0.01*grav*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) WRITE(*,*) 'fq_evapprecip(il,1)=',sigd*0.5*(evap(il,1)+evap(il,2)) WRITE(*,*) 'fq_fluxmasse(il,1)=', 0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il) WRITE(*,*) 'deltaDfq_ddft=',deltaD((xtp(iso_HDO,il,2)-xt(iso_HDO,il,1))/(rp(il,2)-rr(il,1))) WRITE(*,*) 'deltaDfq_evapprecip=',deltaD((xtevap(iso_HDO,il,1)+xtevap(iso_HDO,il,2))/(evap(il,1)+evap(il,2))) WRITE(*,*) 'deltaDfq_fluxmasse=',deltaD((xt(iso_HDO,il,2)-xt(iso_HDO,il,1))/(rr(il,2)-rr(il,1))) WRITE(*,*) 'rr(il,2),rr(il,1)=',rr(il,2),rr(il,1) WRITE(*,*) 'xt(iso_HDO,il,2),xt(iso_HDO,il,1)',xt(iso_HDO,il,2),xt(iso_HDO,il,1) WRITE(*,*) 'dq_tmp=',dq_tmp CALL abort_physic('cv30_routines','cv30_yield',1) endif ! iso_verif_aberrant_enc_nostop endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,1),'cv30_routine 3417') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,1)+delt*fxt(ixt,il,1) enddo IF (iso_verif_tracpos_choix_nostop(xtnew,'cv30_yield 3395',1e-5) & .EQ.1) THEN WRITE(*,*) 'il=',il WRITE(*,*) 'delt,fxt(:,il,1)=',delt,fxt(:,il,1) WRITE(*,*) 'xt(:,il,1)=' ,xt(:,il,1) #ifdef DIAGISO WRITE(*,*) 'fxt_fluxmasse(:,il,1)=',fxt_fluxmasse(:,il,1) WRITE(*,*) 'fxt_ddft(:,il,1)=',fxt_ddft(:,il,1) WRITE(*,*) 'fxt_evapprecip(:,il,1)=', & fxt_evapprecip(:,il,1) WRITE(*,*) 'xt(:,il,2)=',xt(:,il,2) WRITE(*,*) 'xtp(:,il,2)=',xtp(:,il,2) WRITE(*,*) 'xtevap(:,il,1)=',xtevap(:,il,1) WRITE(*,*) 'xtevap(:,il,2)=',xtevap(:,il,2) WRITE(*,*) 'facam,facmp,facev=',0.01*grav*am(il)*work(il), & 0.01*grav*mp(il,2)*work(il),sigd*0.5 #endif ! stop endif #endif #endif ! end cam verif #endif fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il, & 1))+am(il)*(u(il,2)-u(il,1))) fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*(mp(il,2)*(vp(il,2)-v(il, & 1))+am(il)*(v(il,2)-v(il,1))) ELSE ! cvflag_grav fr(il, 1) = 0.1*mp(il, 2)*(rp(il,2)-rr(il,1))*work(il) + & sigd*0.5*(evap(il,1)+evap(il,2)) fr(il, 1) = fr(il, 1) + 0.1*am(il)*(rr(il,2)-rr(il,1))*work(il) #ifdef ISO DO ixt = 1, ntraciso fxt(ixt,il,1)=0.1*mp(il,2)*(xtp(ixt,il,2)-xt(ixt,il,1))*work(il) & +sigd*0.5*(xtevap(ixt,il,1)+xtevap(ixt,il,2)) fxt(ixt,il,1)=fxt(ixt,il,1) & +0.1*am(il)*(xt(ixt,il,2)-xt(ixt,il,1))*work(il) enddo #ifdef DIAGISO fq_ddft(il,1)=fq_ddft(il,1) & +0.1*mp(il,2)*(rp(il,2)-rr(il,1))*work(il) fq_evapprecip(il,1)=fq_evapprecip(il,1) & +sigd*0.5*(evap(il,1)+evap(il,2)) fq_fluxmasse(il,1)=fq_fluxmasse(il,1) & +0.1*am(il)*(rr(il,2)-rr(il,1))*work(il) DO ixt = 1, niso fxt_fluxmasse(ixt,il,1)=fxt(ixt,il,1) & +0.1*am(il)*(xt(ixt,il,2)-xt(ixt,il,1))*work(il) fxt_ddft(ixt,il,1)=fxt(ixt,il,1) & +0.1*mp(il,2)*(xtp(ixt,il,2)-xt(ixt,il,1))*work(il) fxt_evapprecip(ixt,il,1)=fxt(ixt,il,1) & +sigd*0.5*(xtevap(ixt,il,1)+xtevap(ixt,il,2)) enddo #endif ! cam verif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,1), & fr(il,1),'cv30_routines 3023', & errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN IF ((iso_HDO.gt.0).AND. & (rr(il,1)+delt*fr(il,1).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre((xt(iso_HDO,il,1) & +delt*fxt(iso_HDO,il,1)) & /(rr(il,1)+delt*fr(il,1)), & 'cv30_yield 3125b, ddft en 1') endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,1), & 'cv30_routine 3417') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,1)+delt*fxt(ixt,il,1) enddo IF (iso_verif_tracpos_choix_nostop(xtnew, & 'cv30_yield 3449',1e-5) & .EQ.1) THEN WRITE(*,*) 'il=',il WRITE(*,*) 'delt,fxt(:,il,1)=',delt,fxt(:,il,1) WRITE(*,*) 'xt(:,il,1)=' ,xt(:,il,1) ! stop endif #endif #endif ! end cam verif #endif fu(il, 1) = fu(il, 1) + 0.1*work(il)*(mp(il,2)*(up(il,2)-u(il, & 1))+am(il)*(u(il,2)-u(il,1))) fv(il, 1) = fv(il, 1) + 0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il, & 1))+am(il)*(v(il,2)-v(il,1))) END IF ! cvflag_grav END DO ! il ! do j=1,ntra ! do il=1,ncum ! if (cvflag_grav) THEN ! ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il) ! : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) ! : +am(il)*(tra(il,2,j)-tra(il,1,j))) ! else ! ftra(il,1,j)=ftra(il,1,j)+0.1*work(il) ! : *(mp(il,2)*(trap(il,2,j)-tra(il,1,j)) ! : +am(il)*(tra(il,2,j)-tra(il,1,j))) ! END IF ! enddo ! enddo DO j = 2, nl DO il = 1, ncum IF (j<=inb(il)) THEN IF (cvflag_grav) THEN fr(il, 1) = fr(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(qent(il, & j,1)-rr(il,1)) fu(il, 1) = fu(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(uent(il, & j,1)-u(il,1)) fv(il, 1) = fv(il, 1) + 0.01*grav*work(il)*ment(il, j, 1)*(vent(il, & j,1)-v(il,1)) #ifdef ISO DO ixt = 1, ntraciso fxt(ixt,il,1)=fxt(ixt,il,1) & +0.01*grav*work(il)*ment(il,j,1)*(xtent(ixt,il,j,1)-xt(ixt,il,1)) enddo #ifdef DIAGISO fq_detrainement(il,1)=fq_detrainement(il,1) & +0.01*grav*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1)) f_detrainement(il,1)=f_detrainement(il,1) & +0.01*grav*work(il)*ment(il,j,1) q_detrainement(il,1)=q_detrainement(il,1) & +0.01*grav*work(il)*ment(il,j,1)*qent(il,j,1) DO ixt = 1, niso fxt_detrainement(ixt,il,1)=fxt_detrainement(ixt,il,1) & +0.01*grav*work(il)*ment(il,j,1)*(xtent(ixt,il,j,1)-xt(ixt,il,1)) xt_detrainement(ixt,il,1)=xt_detrainement(ixt,il,1) & +0.01*grav*work(il)*ment(il,j,1)*xtent(ixt,il,j,1) enddo #endif ! cam verif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,1), & fr(il,1),'cv30_routines 3251',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN IF ((iso_HDO.gt.0).AND. & (rr(il,1)+delt*fr(il,1).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre((xt(iso_HDO,il,1) & +delt*fxt(iso_HDO,il,1))/(rr(il,1)+delt*fr(il,1)), & 'cv30_yield 3127, dtr melanges') endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,1),'cv30_routine 3417') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,1)+delt*fxt(ixt,il,1) enddo IF (iso_verif_tracpos_choix_nostop(xtnew,'cv30_yield 3525',1e-5) & .EQ.1) THEN WRITE(*,*) 'il=',il WRITE(*,*) 'delt,fxt(:,il,1)=',delt,fxt(:,il,1) WRITE(*,*) 'fac=', 0.01*grav*work(il)*ment(il,j,1) WRITE(*,*) 'xt(:,il,1)=' ,xt(:,il,1) WRITE(*,*) 'xtent(:,il,j,1)=' ,xtent(:,il,j,1) ! stop endif #endif #endif ! end cam verif #endif ELSE ! cvflag_grav fr(il, 1) = fr(il, 1) + 0.1*work(il)*ment(il, j, 1)*(qent(il,j,1)- & rr(il,1)) fu(il, 1) = fu(il, 1) + 0.1*work(il)*ment(il, j, 1)*(uent(il,j,1)-u & (il,1)) fv(il, 1) = fv(il, 1) + 0.1*work(il)*ment(il, j, 1)*(vent(il,j,1)-v & (il,1)) #ifdef ISO DO ixt = 1, ntraciso fxt(ixt,il,1)=fxt(ixt,il,1) & +0.1*work(il)*ment(il,j,1)*(xtent(ixt,il,j,1)-xt(ixt,il,1)) enddo #ifdef DIAGISO fq_detrainement(il,1)=fq_detrainement(il,1) & +0.1*work(il)*ment(il,j,1)*(qent(il,j,1)-rr(il,1)) f_detrainement(il,1)=f_detrainement(il,1) & +0.1*work(il)*ment(il,j,1) q_detrainement(il,1)=q_detrainement(il,1) & +0.1*work(il)*ment(il,j,1)*qent(il,j,1) DO ixt = 1, niso fxt_detrainement(ixt,il,1)=fxt_detrainement(ixt,il,1) & +0.1*work(il)*ment(il,j,1)*(xtent(ixt,il,j,1)-xt(ixt,il,1)) xt_detrainement(ixt,il,1)=xt_detrainement(ixt,il,1) & +0.1*work(il)*ment(il,j,1)*xtent(ixt,il,j,1) enddo #endif ! cam verif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,1), & fr(il,1),'cv30_routines 3092',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN IF ((iso_HDO.gt.0).AND. & (rr(il,1)+delt*fr(il,1).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre((xt(iso_HDO,il,1) & +delt*fxt(iso_HDO,il,1))/(rr(il,1)+delt*fr(il,1)), & 'cv30_yield 3127b, dtr melanges') endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,1),'cv30_routine 3462') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,1)+delt*fxt(ixt,il,1) enddo IF (iso_verif_tracpos_choix_nostop(xtnew,'cv30_yield 3753',1e-5) & .EQ.1) THEN WRITE(*,*) 'il=',il endif #endif #endif ! end cam verif #endif END IF ! cvflag_grav END IF ! j END DO END DO ! do k=1,ntra ! do j=2,nl ! do il=1,ncum ! if (j.le.inb(il)) THEN ! if (cvflag_grav) THEN ! ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1) ! : *(traent(il,j,1,k)-tra(il,1,k)) ! else ! ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1) ! : *(traent(il,j,1,k)-tra(il,1,k)) ! END IF ! END IF ! enddo ! enddo ! enddo ! *** calculate tendencies of potential temperature and mixing ratio *** ! *** at levels above the lowest level *** ! *** first find the net saturated updraft and downdraft mass fluxes *** ! *** through each level *** DO i = 2, nl + 1 ! newvecto: mettre nl au lieu nl+1? num1 = 0 DO il = 1, ncum IF (i<=inb(il)) num1 = num1 + 1 END DO IF (num1<=0) GO TO 500 CALL zilch(amp1, ncum) CALL zilch(ad, ncum) DO k = i + 1, nl + 1 DO il = 1, ncum IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN amp1(il) = amp1(il) + m(il, k) END IF END DO END DO DO k = 1, i DO j = i + 1, nl + 1 DO il = 1, ncum IF (i<=inb(il) .AND. j<=(inb(il)+1)) THEN amp1(il) = amp1(il) + ment(il, k, j) END IF END DO END DO END DO DO k = 1, i - 1 DO j = i, nl + 1 ! newvecto: nl au lieu nl+1? DO il = 1, ncum IF (i<=inb(il) .AND. j<=inb(il)) THEN ad(il) = ad(il) + ment(il, j, k) END IF END DO END DO END DO DO il = 1, ncum IF (i<=inb(il)) THEN dpinv = 1.0/(ph(il,i)-ph(il,i+1)) cpinv = 1.0/cpn(il, i) ! convect3 if((0.1*dpinv*amp1).ge.delti)iflag(il)=4 IF (cvflag_grav) THEN IF ((0.01*grav*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto ELSE IF ((0.1*dpinv*amp1(il))>=delti) iflag(il) = 1 ! vecto END IF IF (cvflag_grav) THEN ft(il, i) = 0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il, & i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, & i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) - 0.5*sigd*lvcp(il, i)*(evap( & il,i)+evap(il,i+1)) rat = cpn(il, i-1)*cpinv ft(il, i) = ft(il, i) - 0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i) & -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv ft(il, i) = ft(il, i) + 0.01*grav*dpinv*ment(il, i, i)*(hp(il,i)-h( & il,i)+t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv ELSE ! cvflag_grav ft(il, i) = 0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il, & i)+(gz(il,i+1)-gz(il,i))*cpinv)-ad(il)*(t(il,i)-t(il, & i-1)+(gz(il,i)-gz(il,i-1))*cpinv)) - 0.5*sigd*lvcp(il, i)*(evap( & il,i)+evap(il,i+1)) rat = cpn(il, i-1)*cpinv ft(il, i) = ft(il, i) - 0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i)-mp(il & ,i)*t(il,i-1)*rat*b(il,i-1))*dpinv ft(il, i) = ft(il, i) + 0.1*dpinv*ment(il, i, i)*(hp(il,i)-h(il,i)+ & t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv END IF ! cvflag_grav ft(il, i) = ft(il, i) + 0.01*sigd*wt(il, i)*(cl-cpd)*water(il, i+1)*( & t(il,i+1)-t(il,i))*dpinv*cpinv IF (cvflag_grav) THEN fr(il, i) = 0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il, & i))-ad(il)*(rr(il,i)-rr(il,i-1))) fu(il, i) = fu(il, i) + 0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il, & i))-ad(il)*(u(il,i)-u(il,i-1))) fv(il, i) = fv(il, i) + 0.01*grav*dpinv*(amp1(il)*(v(il,i+1)-v(il, & i))-ad(il)*(v(il,i)-v(il,i-1))) #ifdef ISO #ifdef DIAGISO fq_fluxmasse(il,i)=fq_fluxmasse(il,i) & +0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) & -ad(il)*(rr(il,i)-rr(il,i-1))) ! modif 2 fev: pour avoir subsidence compensatoire totale, on retranche ! ad. #endif ! ici, on separe 2 cas, pour eviter le cas pathologique decrit plus haut ! pour la tendance liee a Am en i=1, qui peut conduire a des isotopes ! negatifs dans les cas ou les flux de masse soustrait plus de 90% de la ! vapeur de la couche. Voir plus haut le detail des equations. ! La difference ici est qu'on considere les flux de masse amp1 et ad en ! meme temps. dq_tmp= 0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) & -ad(il)*(rr(il,i)-rr(il,i-1)))*delt ! c'est equivalent a dqi= kamp1*qip1+kad*qim1-(kamp1+kad)*qi IF ((dq_tmp/rr(il,i).lt.-0.9).AND.correction_excess_aberrant) THEN ! nouvelle formulation k_tmp=0.01*grav*dpinv*amp1(il)*delt kad_tmp=0.01*grav*dpinv*ad(il)*delt DO ixt = 1, ntraciso R_tmp=(xt(ixt,il,i)+k_tmp*xt(ixt,il,i+1)+kad_tmp*xt(ixt,il,i-1)) & /(rr(il,i)+k_tmp*rr(il,i+1)+kad_tmp*rr(il,i-1)) dx_tmp= R_tmp*( rr(il,i)+ dq_tmp)-xt(ixt,il,i) fxt(ixt,il,i)= dx_tmp/delt #ifdef ISOVERIF IF ((ixt.EQ.iso_HDO).OR.(ixt.EQ.iso_eau)) THEN WRITE(*,*) 'cv30_routines 4367: il,i,ixt=',il,i,ixt WRITE(*,*) 'dq_tmp,rr(il,i)=',dq_tmp,rr(il,i) WRITE(*,*) 'R_tmp,dx_tmp,delt=',R_tmp,dx_tmp,delt WRITE(*,*) 'amp1(il),ad(il)=',amp1(il),ad(il) WRITE(*,*) 'xt(ixt,il,i-1:i+1)=',xt(ixt,il,i-1:i+1) WRITE(*,*) 'rr(il,i-1:i+1)=',rr(il,i-1:i+1) WRITE(*,*) 'fxt=',dx_tmp/delt WRITE(*,*) 'rr(il,i)+dq_tmp=',rr(il,1)+dq_tmp WRITE(*,*) 'xt(ixt,il,i)+dx_tmp=',xt(ixt,il,i)+dx_tmp WRITE(*,*) 'xt(ixt,il,i)+fxt(ixt,il,i)*delt=', & xt(ixt,il,i)+fxt(ixt,il,i)*delt WRITE(*,*) 'fxt(:,il,i)=',fxt(:,il,i) endif !if (ixt.EQ.iso_HDO) THEN #endif enddo ! do ixt = 1, ntraciso #ifdef DIAGISO DO ixt = 1, niso fxt_fluxmasse(ixt,il,i)=fxt(ixt,il,i) enddo #endif else !if (dq_tmp/rr(il,i).lt.-0.9) THEN ! ancienne formulation DO ixt = 1, ntraciso fxt(ixt,il,i)= & 0.01*grav*dpinv*(amp1(il)*(xt(ixt,il,i+1)-xt(ixt,il,i)) & -ad(il)*(xt(ixt,il,i)-xt(ixt,il,i-1))) enddo #ifdef DIAGISO DO ixt = 1, niso fxt_fluxmasse(ixt,il,i)=fxt_fluxmasse(ixt,il,i)+ & 0.01*grav*dpinv*(amp1(il)*(xt(ixt,il,i+1)-xt(ixt,il,i)) & -ad(il)*(xt(ixt,il,i)-xt(ixt,il,i-1))) enddo #endif endif !if (dq_tmp/rr(il,i).lt.-0.9) THEN ! cam verif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,i), & fr(il,i),'cv30_routines 3226',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN DO ixt=1,niso CALL iso_verif_noNAN(fxt(ixt,il,i),'cv30_routines 3229') enddo IF ((iso_HDO.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre((xt(iso_HDO,il,i) & +delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 3384, flux masse') endif !if (iso_HDO.gt.0) THEN IF ((iso_HDO.gt.0).AND.(iso_O18.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*fxt(iso_HDO,il,i))/(rr(il,i)+delt*fr(il,i)), & (xt(iso_O18,il,i)+delt*fxt(iso_O18,il,i))/(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 3384,O18, flux masse') endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,1),'cv30_routine 3626') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,i)+delt*fxt(ixt,il,i) enddo IF (iso_verif_tracpos_choix_nostop(xtnew,'cv30_yield 3727',1e-5) & .EQ.1) THEN WRITE(*,*) 'il,i=',il,i WRITE(*,*) 'fxt(:,il,i)=',fxt(:,il,i) WRITE(*,*) 'amp1(il),ad(il),fac=', & amp1(il),ad(il),0.01*grav*dpinv WRITE(*,*) 'xt(:,il,i+1)=' ,xt(:,il,i+1) WRITE(*,*) 'xt(:,il,i)=' ,xt(:,il,i) WRITE(*,*) 'xt(:,il,i-1)=' ,xt(:,il,i-1) ! stop endif #endif #endif ! end cam verif #endif ELSE ! cvflag_grav fr(il, i) = 0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il, & i))-ad(il)*(rr(il,i)-rr(il,i-1))) fu(il, i) = fu(il, i) + 0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il, & i))-ad(il)*(u(il,i)-u(il,i-1))) fv(il, i) = fv(il, i) + 0.1*dpinv*(amp1(il)*(v(il,i+1)-v(il, & i))-ad(il)*(v(il,i)-v(il,i-1))) #ifdef ISO DO ixt = 1, ntraciso fxt(ixt,il,i)= & 0.1*dpinv*(amp1(il)*(xt(ixt,il,i+1)-xt(ixt,il,i)) & -ad(il)*(xt(ixt,il,i)-xt(ixt,il,i-1))) enddo #ifdef DIAGISO fq_fluxmasse(il,i)=fq_fluxmasse(il,i) & +0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i)) & -ad(il)*(rr(il,i)-rr(il,i-1))) DO ixt = 1, niso fxt_fluxmasse(ixt,il,i)=fxt_fluxmasse(ixt,il,i)+ & 0.1*dpinv*(amp1(il)*(xt(ixt,il,i+1)-xt(ixt,il,i)) & -ad(il)*(xt(ixt,il,i)-xt(ixt,il,i-1))) enddo #endif ! cam verif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,i), & fr(il,i),'cv30_routines 3252',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN DO ixt=1,niso CALL iso_verif_noNAN(fxt(ixt,il,i),'cv30_routines 3229') enddo ! correction 21 oct 2008 IF ((iso_HDO.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre((xt(iso_HDO,il,i) & +delt*fxt(iso_HDO,il,i))/(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 3384b flux masse') IF (iso_O18.gt.0) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)), & (xt(iso_O18,il,i)+delt*fxt(iso_O18,il,i)) & /(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 3384bO18 flux masse') endif !if (iso_O18.gt.0) THEN endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,1),'cv30_routine 3674') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,i)+delt*fxt(ixt,il,i) enddo IF (iso_verif_tracpos_choix_nostop(xtnew,'cv30_yield 3775',1e-5) & .EQ.1) THEN WRITE(*,*) 'il,i=',il,i endif #endif #endif ! end cam verif #endif END IF ! cvflag_grav END IF ! i END DO ! do k=1,ntra ! do il=1,ncum ! if (i.le.inb(il)) THEN ! dpinv=1.0/(ph(il,i)-ph(il,i+1)) ! cpinv=1.0/cpn(il,i) ! if (cvflag_grav) THEN ! ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv ! : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) ! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) ! else ! ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv ! : *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k)) ! : -ad(il)*(tra(il,i,k)-tra(il,i-1,k))) ! END IF ! END IF ! enddo ! enddo DO k = 1, i - 1 DO il = 1, ncum IF (i<=inb(il)) THEN dpinv = 1.0/(ph(il,i)-ph(il,i+1)) cpinv = 1.0/cpn(il, i) awat = elij(il, k, i) - (1.-ep(il,i))*clw(il, i) awat = amax1(awat, 0.0) #ifdef ISO ! on change le traitement de cette ligne le 8 mai 2009: ! avant, on avait: xtawat=xtelij(il,k,i)-(1.-xtep(il,i))*xtclw(il,i) ! c'est a dire que Rawat=Relij+(1-ep)*clw/awat*(Relij-Rclw) ! si Relij!=Rclw, alors un fractionnement isotopique non physique etait ! introduit. ! En fait, awat represente le surplus de condensat dans le melange par ! rapport a celui restant dans la colonne adiabatique ! ce surplus a la meme compo que le elij, sans fractionnement. ! d'ou le nouveau traitement ci-dessous. IF (elij(il,k,i).gt.0.0) THEN DO ixt = 1, ntraciso xtawat(ixt)=awat*(xtelij(ixt,il,k,i)/elij(il,k,i)) ! xtawat(ixt)=amax1(xtawat(ixt),0.0) ! pas necessaire enddo else !if (elij(il,k,i).gt.0.0) THEN ! normalement, si elij(il,k,i)<=0, alors awat=0 ! on le verifie. Si c'est vrai -> xtawat=0 aussi #ifdef ISOVERIF CALL iso_verif_egalite(awat,0.0,'cv30_yield 3779') #endif DO ixt = 1, ntraciso xtawat(ixt)=0.0 enddo endif ! cam verif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(xtawat(iso_eau), & awat,'cv30_routines 3301',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(xtawat(1),'cv30_routine 3729') #endif #endif ! end cam verif #endif IF (cvflag_grav) THEN fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k & ,i)-awat-rr(il,i)) fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k & ,i)-u(il,i)) fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k & ,i)-v(il,i)) #ifdef ISO DO ixt = 1, ntraciso fxt(ixt,il,i)=fxt(ixt,il,i) & +0.01*grav*dpinv*ment(il,k,i) & *(xtent(ixt,il,k,i)-xtawat(ixt)-xt(ixt,il,i)) enddo #ifdef DIAGISO fq_detrainement(il,i)=fq_detrainement(il,i) & +0.01*grav*dpinv*ment(il,k,i) & *(qent(il,k,i)-awat-rr(il,i)) f_detrainement(il,i)=f_detrainement(il,i)& +0.01*grav*dpinv*ment(il,k,i) q_detrainement(il,i)=q_detrainement(il,i) & +0.01*grav*dpinv*ment(il,k,i)*qent(il,k,i) DO ixt = 1, niso fxt_detrainement(ixt,il,i)=fxt_detrainement(ixt,il,i) & +0.01*grav*dpinv*ment(il,k,i) & *(xtent(ixt,il,k,i)-xtawat(ixt)-xt(ixt,il,i)) xt_detrainement(ixt,il,i)=xt_detrainement(ixt,il,i) & +0.01*grav*dpinv*ment(il,k,i)*xtent(ixt,il,k,i) enddo #endif ! cam verif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,i), & fr(il,i),'cv30_routines 3325',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN DO ixt=1,niso CALL iso_verif_noNAN(fxt(ixt,il,i),'cv30_routines 3328') enddo IF ((iso_HDO.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN IF (iso_verif_aberrant_enc_nostop((xt(iso_HDO,il,i) & +delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)),'cv30_yield 3396a, dtr mels') & .EQ.1) THEN WRITE(*,*) 'il,k,i=',il,k,i WRITE(*,*) 'rr,delt,fr=',rr(il,i),delt,fr(il,i) WRITE(*,*) 'frnew=',0.01*grav*dpinv*ment(il, k, i)*(qent(il,k,i)-awat-rr(il,i)) WRITE(*,*) 'frold=',fr(il,i)-0.01*grav*dpinv*ment(il, k, i)*(qent(il,k,i)-awat-rr(il,i)) WRITE(*,*) 'deltaDfrnew=',deltaD((xtent(iso_HDO,il,k,i)-xtawat(iso_HDO)-xt(iso_HDO,il,i)) & /(qent(il,k,i)-awat-rr(il,i))) WRITE(*,*) 'deltaDfrold=',deltaD((fxt(iso_HDO,il,i) & -0.01*grav*dpinv*ment(il, k, i)*(xtent(iso_HDO,il,k,i)-xtawat(iso_HDO)-xt(iso_HDO,il,i))) & /(fr(il,i)-0.01*grav*dpinv*ment(il, k, i)*(qent(il,k,i)-awat-rr(il,i)))) WRITE(*,*) 'q+=',rr(il,i)+delt*fr(il,i) WRITE(*,*) 'qent,awat=',qent(il,k,i),awat WRITE(*,*) 'elij,clw,ep=',elij(il,k,i),clw(il,i),ep(il,i) WRITE(*,*) 'deltaDfr=',deltaD(fxt(iso_hdo,il,i)/fr(il,i)) WRITE(*,*) 'deltaDrr=',deltaD(xt(iso_hdo,il,i)/rr(il,i)) WRITE(*,*) 'deltaDqent=',deltaD(xtent(iso_hdo,il,k,i) & /qent(il,k,i)) WRITE(*,*) 'deltaDqent-awat=',deltaD((xtent(iso_hdo,il,k,i)-xtawat(iso_HDO)) & /(qent(il,k,i)-awat)) WRITE(*,*) 'deltaDawat=',deltaD(xtawat(iso_hdo)/awat) WRITE(*,*) 'deltaDclw=',deltaD(xtclw(iso_hdo,il,i)/clw(il,i)) ! stop endif IF (iso_O18.gt.0) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)), & (xt(iso_O18,il,i)+delt*fxt(iso_O18,il,i)) & /(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 3396aO18, dtr mels') endif !if (iso_O18.gt.0) THEN endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,i),'cv30_routine 3784') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,i)+delt*fxt(ixt,il,i) enddo IF (iso_verif_tracpos_choix_nostop(xtnew,'cv30_yield 3905',1e-5) & .EQ.1) THEN WRITE(*,*) 'il,i=',il,i endif ! CALL iso_verif_tracpos_choix(xtnew,'cv30_yield 3905',1e-5) #endif #endif #endif ELSE ! cvflag_grav fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)- & awat-rr(il,i)) fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k & ,i)-u(il,i)) fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( & il,i)) #ifdef ISO DO ixt = 1, ntraciso fxt(ixt,il,i)=fxt(ixt,il,i) & +0.1*dpinv*ment(il,k,i) & *(xtent(ixt,il,k,i)-xtawat(ixt)-xt(ixt,il,i)) enddo #ifdef DIAGISO fq_detrainement(il,i)=fq_detrainement(il,i) & +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i)) f_detrainement(il,i)=f_detrainement(il,i) & +0.1*dpinv*ment(il,k,i) q_detrainement(il,i)=q_detrainement(il,i) & +0.1*dpinv*ment(il,k,i)*qent(il,k,i) DO ixt = 1, niso fxt_detrainement(ixt,il,i)=fxt_detrainement(ixt,il,i) & +0.1*dpinv*ment(il,k,i) & *(xtent(ixt,il,k,i)-xtawat(ixt)-xt(ixt,il,i)) xt_detrainement(ixt,il,i)=xt_detrainement(ixt,il,i) & +0.1*dpinv*ment(il,k,i)*xtent(ixt,il,k,i) enddo #endif ! cam verif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,i), & fr(il,i),'cv30_routines 3350',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN DO ixt=1,niso CALL iso_verif_noNAN(fxt(ixt,il,i),'cv30_routines 3353') enddo IF ((iso_HDO.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre((xt(iso_HDO,il,i) & +delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)),'cv30_yield 3396b, dtr mels') endif !if (iso_HDO.gt.0) THEN IF ((iso_HDO.gt.0).AND.(iso_O18.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*fxt(iso_HDO,il,i))/(rr(il,i)+delt*fr(il,i)), & (xt(iso_O18,il,i)+delt*fxt(iso_O18,il,i))/(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 3396b,O18, dtr mels') endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,i),'cv30_routine 3828') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,i)+delt*fxt(ixt,il,i) enddo IF (iso_verif_tracpos_choix_nostop(xtnew,'cv30_yield 3949',1e-5) & .EQ.1) THEN WRITE(*,*) 'il,i=',il,i endif ! CALL iso_verif_tracpos_choix(xtnew,'cv30_yield 3949',1e-5) #endif #endif ! end cam verif #endif END IF ! cvflag_grav ! (saturated updrafts resulting from mixing) ! cld qcond(il, i) = qcond(il, i) + (elij(il,k,i)-awat) ! cld nqcond(il, i) = nqcond(il, i) + 1. ! cld END IF ! i END DO END DO ! do j=1,ntra ! do k=1,i-1 ! do il=1,ncum ! if (i.le.inb(il)) THEN ! dpinv=1.0/(ph(il,i)-ph(il,i+1)) ! cpinv=1.0/cpn(il,i) ! if (cvflag_grav) THEN ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) ! : *(traent(il,k,i,j)-tra(il,i,j)) ! else ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) ! : *(traent(il,k,i,j)-tra(il,i,j)) ! END IF ! END IF ! enddo ! enddo ! enddo DO k = i, nl + 1 DO il = 1, ncum IF (i<=inb(il) .AND. k<=inb(il)) THEN dpinv = 1.0/(ph(il,i)-ph(il,i+1)) cpinv = 1.0/cpn(il, i) IF (cvflag_grav) THEN fr(il, i) = fr(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(qent(il,k & ,i)-rr(il,i)) fu(il, i) = fu(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(uent(il,k & ,i)-u(il,i)) fv(il, i) = fv(il, i) + 0.01*grav*dpinv*ment(il, k, i)*(vent(il,k & ,i)-v(il,i)) #ifdef ISO DO ixt = 1, ntraciso fxt(ixt,il,i)=fxt(ixt,il,i) & +0.01*grav*dpinv*ment(il,k,i)*(xtent(ixt,il,k,i)-xt(ixt,il,i)) enddo #ifdef DIAGISO fq_detrainement(il,i)=fq_detrainement(il,i) & +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i)) f_detrainement(il,i)=f_detrainement(il,i) & +0.01*grav*dpinv*ment(il,k,i) q_detrainement(il,i)=q_detrainement(il,i) & +0.01*grav*dpinv*ment(il,k,i)*qent(il,k,i) DO ixt = 1, niso fxt_detrainement(ixt,il,i)=fxt_detrainement(ixt,il,i) & +0.01*grav*dpinv*ment(il,k,i)*(xtent(ixt,il,k,i)-xt(ixt,il,i)) xt_detrainement(ixt,il,i)=xt_detrainement(ixt,il,i) & +0.01*grav*dpinv*ment(il,k,i)*xtent(ixt,il,k,i) enddo #endif ! cam verif #ifdef ISOVERIF IF ((il.EQ.1636).AND.(i.EQ.9)) THEN WRITE(*,*) 'cv30 4785: on ajoute le dtr ici:' WRITE(*,*) 'M=',0.01*grav*dpinv*ment(il, k, i) WRITE(*,*) 'q,qe=',rr(il,i),qent(il,k,i) bx=0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i)) DO ixt=1,niso xtbx(ixt)=0.01*grav*dpinv*ment(il,k,i)*(xtent(ixt,il,k,i)-xt(ixt,il,i)) enddo endif DO ixt=1,niso CALL iso_verif_noNaN(fxt(ixt,il,i),'cv30_yield 4351') enddo #endif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,i), & fr(il,i),'cv30_routines 3408',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN DO ixt=1,niso CALL iso_verif_noNAN(fxt(ixt,il,i),'cv30_routines 3411') enddo IF (1.EQ.0) THEN IF ((iso_HDO.gt.0).AND.(delt*fr(il,i).gt.ridicule)) THEN IF (iso_verif_aberrant_enc_nostop( & fxt(iso_HDO,il,i)/fr(il,i), & 'cv30_yield 3572, dtr mels').EQ.1) THEN WRITE(*,*) 'i,icb(il),inb(il)=',i,icb(il),inb(il) WRITE(*,*) 'fr(il,i)=',fr(il,i) ! if (fr(il,i).gt.ridicule*1e5) THEN ! stop ! endif endif endif !if (iso_HDO.gt.0) THEN endif !if (1.EQ.0) THEN IF ((iso_HDO.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre((xt(iso_HDO,il,i) & +delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)),'cv30_yield 3605, dtr mels') IF (iso_O18.gt.0) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)), & (xt(iso_O18,il,i)+delt*fxt(iso_O18,il,i)) & /(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 3605O18, dtr mels') IF ((il.EQ.1636).AND.(i.EQ.9)) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*(fxt(iso_HDO,il,i)-xtbx(iso_HDO))) & /(rr(il,i)+delt*(fr(il,i)-bx)), & (xt(iso_O18,il,i)+delt*(fxt(iso_O18,il,i)-xtbx(iso_O18))) & /(rr(il,i)+delt*(fr(il,i)-bx)), & 'cv30_yield 3605O18_nobx, dtr mels') endif !if ((il.EQ.1636).AND.(i.EQ.9)) THEN endif !if (iso_O18.gt.0) THEN endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,i),'cv30_routine 3921') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,i)+delt*fxt(ixt,il,i) enddo IF (iso_verif_tracpos_choix_nostop(xtnew,'cv30_yield 4036',1e-5) & .EQ.1) THEN WRITE(*,*) 'il,i=',il,i endif ! CALL iso_verif_tracpos_choix(xtnew,'cv30_yield 4036',1e-5) #endif #endif ! end cam verif #endif ELSE ! cvflag_grav fr(il, i) = fr(il, i) + 0.1*dpinv*ment(il, k, i)*(qent(il,k,i)-rr & (il,i)) fu(il, i) = fu(il, i) + 0.1*dpinv*ment(il, k, i)*(uent(il,k,i)-u( & il,i)) fv(il, i) = fv(il, i) + 0.1*dpinv*ment(il, k, i)*(vent(il,k,i)-v( & il,i)) #ifdef ISO DO ixt = 1, ntraciso fxt(ixt,il,i)=fxt(ixt,il,i) & +0.1*dpinv*ment(il,k,i)*(xtent(ixt,il,k,i)-xt(ixt,il,i)) enddo #ifdef DIAGISO fq_detrainement(il,i)=fq_detrainement(il,i) & +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i)) f_detrainement(il,i)=f_detrainement(il,i) & +0.1*dpinv*ment(il,k,i) q_detrainement(il,i)=q_detrainement(il,i) & +0.1*dpinv*ment(il,k,i)*qent(il,k,i) DO ixt = 1, niso fxt_detrainement(ixt,il,i)=fxt_detrainement(ixt,il,i) & +0.1*dpinv*ment(il,k,i)*(xtent(ixt,il,k,i)-xt(ixt,il,i)) xt_detrainement(ixt,il,i)=xt_detrainement(ixt,il,i) & +0.1*dpinv*ment(il,k,i)*xtent(ixt,il,k,i) enddo #endif ! cam verif #ifdef ISOVERIF IF ((il.EQ.1636).AND.(i.EQ.9)) THEN WRITE(*,*) 'cv30 4785b: on ajoute le dtr ici:' WRITE(*,*) 'M=',0.1*dpinv*ment(il, k, i) WRITE(*,*) 'q,qe=',rr(il,i),qent(il,k,i) endif IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,i), & fr(il,i),'cv30_routines 3433',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN DO ixt=1,niso CALL iso_verif_noNAN(fxt(ixt,il,i),'cv30_routines 3436') enddo IF ((iso_HDO.gt.0).AND.(delt*fr(il,i).gt.ridicule)) THEN IF (iso_verif_aberrant_enc_nostop( & fxt(iso_HDO,il,i)/fr(il,i), & 'cv30_yield 3597').EQ.1) THEN WRITE(*,*) 'i,icb(il),inb(il)=',i,icb(il),inb(il) stop endif endif !if (iso_HDO.gt.0) THEN IF ((iso_HDO.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre((xt(iso_HDO,il,i) & +delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)),'cv30_yield 3605b, dtr mels') endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,i),'cv30_routine 3972') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,i)+delt*fxt(ixt,il,i) enddo IF (iso_verif_tracpos_choix_nostop(xtnew,'cv30_yield 4091',1e-5) & .EQ.1) THEN WRITE(*,*) 'il,i=',il,i endif ! CALL iso_verif_tracpos_choix(xtnew,'cv30_yield 4091',1e-5) #endif #endif ! end cam verif #endif END IF ! cvflag_grav END IF ! i and k END DO END DO ! do j=1,ntra ! do k=i,nl+1 ! do il=1,ncum ! if (i.le.inb(il) .AND. k.le.inb(il)) THEN ! dpinv=1.0/(ph(il,i)-ph(il,i+1)) ! cpinv=1.0/cpn(il,i) ! if (cvflag_grav) THEN ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i) ! : *(traent(il,k,i,j)-tra(il,i,j)) ! else ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i) ! : *(traent(il,k,i,j)-tra(il,i,j)) ! END IF ! END IF ! i and k ! enddo ! enddo ! enddo DO il = 1, ncum IF (i<=inb(il)) THEN dpinv = 1.0/(ph(il,i)-ph(il,i+1)) cpinv = 1.0/cpn(il, i) IF (cvflag_grav) THEN ! sb: on ne fait pas encore la correction permettant de mieux ! conserver l'eau: fr(il, i) = fr(il, i) + 0.5*sigd*(evap(il,i)+evap(il,i+1)) + & 0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)*(rp(il, & i)-rr(il,i-1)))*dpinv fu(il, i) = fu(il, i) + 0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il, & i))-mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv fv(il, i) = fv(il, i) + 0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il, & i))-mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv #ifdef ISO DO ixt = 1, niso fxt(ixt,il,i)=fxt(ixt,il,i) & +0.5*sigd*(xtevap(ixt,il,i)+xtevap(ixt,il,i+1)) & +0.01*grav*(mp(il,i+1)*(xtp(ixt,il,i+1)-xt(ixt,il,i)) & -mp(il,i) & *(xtp(ixt,il,i)-xt(ixt,il,i-1)))*dpinv enddo #ifdef DIAGISO fq_evapprecip(il,i)=fq_evapprecip(il,i) & +0.5*sigd*(evap(il,i)+evap(il,i+1)) fq_ddft(il,i)=fq_ddft(il,i) & +0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) & *(rp(il,i)-rr(il,i-1)))*dpinv DO ixt = 1, niso fxt_evapprecip(ixt,il,i)=fxt_evapprecip(ixt,il,i) & +0.5*sigd*(xtevap(ixt,il,i)+xtevap(ixt,il,i+1)) fxt_ddft(ixt,il,i)=fxt_ddft(ixt,il,i) & +0.01*grav*(mp(il,i+1)*(xtp(ixt,il,i+1)-xt(ixt,il,i)) & -mp(il,i)*(xtp(ixt,il,i)-xt(ixt,il,i-1)))*dpinv enddo #endif #ifdef ISOVERIF DO ixt=1,niso CALL iso_verif_noNaN(xt(ixt,il,i),'cv30_yield 4514') CALL iso_verif_noNaN(fxt(ixt,il,i),'cv30_yield 4515') enddo IF ((iso_HDO.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN IF (iso_verif_aberrant_enc_nostop((xt(iso_HDO,il,i) & +delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)),'cv30_yield 4175') & .EQ.1) THEN WRITE(*,*) 'il,i=',il,i IF (rr(il,i).NE.0.0) THEN WRITE(*,*) 'il,i,rr,deltaD=',il,i,rr(il,i),deltaD & (xt(iso_HDO,il,i)/rr(il,i)) endif IF (fr(il,i).NE.0.0) THEN WRITE(*,*) 'fr,fxt,deltaD=',fr(il,i),fxt(iso_HDO,il,i), & deltaD(fxt(iso_HDO,il,i)/fr(il,i)) endif #ifdef DIAGISO IF (fq_ddft(il,i).NE.0.0) THEN WRITE(*,*) 'fq_ddft,deltaD=',fq_ddft(il,i),deltaD( & fxt_ddft(iso_HDO,il,i)/fq_ddft(il,i)) endif IF (fq_evapprecip(il,i).NE.0.0) THEN WRITE(*,*) 'fq_evapprecip,deltaD=',fq_evapprecip(il,i),deltaD( & fxt_evapprecip(iso_HDO,il,i)/fq_evapprecip(il,i)) endif #endif WRITE(*,*) 'sigd,evap(il,i),evap(il,i+1)=', & sigd,evap(il,i),evap(il,i+1) WRITE(*,*) 'xtevap(iso_HDO,il,i),xtevap(iso_HDO,il,i+1)=', & xtevap(iso_HDO,il,i),xtevap(iso_HDO,il,i+1) WRITE(*,*) 'grav,mp(il,i+1),mp(il,i),dpinv=', & grav,mp(il,i+1),mp(il,i),dpinv WRITE(*,*) 'rp(il,i+1),rr(il,i),rp(il,i),rr(il,i-1)=', & rp(il,i+1),rr(il,i),rp(il,i),rr(il,i-1) WRITE(*,*) 'xtp(il,i+1),xt(il,i),xtp(il,i),xt(il,i-1)=', & xtp(iso_HDO,il,i+1),xt(iso_HDO,il,i), & xtp(iso_HDO,il,i),xt(iso_HDO,il,i-1) stop endif endif !if (iso_HDO.gt.0) THEN IF ((iso_HDO.gt.0).AND.(iso_O18.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*fxt(iso_HDO,il,i))/(rr(il,i)+delt*fr(il,i)), & (xt(iso_O18,il,i)+delt*fxt(iso_O18,il,i))/(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 5029,O18, evap') IF ((il.EQ.1636).AND.(i.EQ.9)) THEN WRITE(*,*) 'cv30_yield 5057: ici, on verifie deltaD_nobx' WRITE(*,*) 'il,i=',il,i WRITE(*,*) 'fr(il,i),bx,fr(il,i)-bx=',fr(il,i),bx,fr(il,i)-bx WRITE(*,*) 'q,q+=',rr(il,i),rr(il,i)+delt*(fr(il,i)-bx) WRITE(*,*) 'deltaD,deltaD+=',deltaD(xt(iso_HDO,il,inb(il))/rr(il,inb(il))), & deltaD( (xt(iso_HDO,il,i)+delt*(fxt(iso_HDO,il,i)-xtbx(iso_HDO)))/(rr(il,i)+delt*(fr(il,i)-bx))) WRITE(*,*) 'deltaO18,deltaO18+=',deltaO(xt(iso_O18,il,inb(il))/rr(il,inb(il))), & deltaO( (xt(iso_O18,il,i)+delt*(fxt(iso_O18,il,i)-xtbx(iso_O18)))/(rr(il,i)+delt*(fr(il,i)-bx))) CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*(fxt(iso_HDO,il,i)-xtbx(iso_HDO))) & /(rr(il,i)+delt*(fr(il,i)-bx)), & (xt(iso_O18,il,i)+delt*(fxt(iso_O18,il,i)-xtbx(iso_O18))) & /(rr(il,i)+delt*(fr(il,i)-bx)), & 'cv30_yield 5029_nobx,O18, evap, no bx') endif !if ((il.EQ.1636).AND.(i.EQ.9)) THEN endif !if (iso_HDO.gt.0) THEN #endif #ifdef ISOTRAC IF ((option_traceurs.NE.6).AND.(option_traceurs.NE.19)) THEN ! facile: on fait comme l'eau DO ixt = 1+niso,ntraciso fxt(ixt,il,i)=fxt(ixt,il,i) & +0.5*sigd*(xtevap(ixt,il,i)+xtevap(ixt,il,i+1)) & +0.01*grav*(mp(il,i+1)*(xtp(ixt,il,i+1)-xt(ixt,il,i)) & -mp(il,i)*(xtp(ixt,il,i)-xt(ixt,il,i-1)))*dpinv enddo !do ixt = 1+niso,ntraciso else ! taggage des ddfts: ! la formule pour fq_ddft suppose que le ddft est en RP. Ce n'est pas le ! cas pour le water tagging puisqu'il y a conversion des molecules ! blances entrainees en molecule rouges. ! Il faut donc prendre en compte ce taux de conversion quand ! entrainement d'env vers ddft ! conversion(iiso)=0.01*grav*dpinv ! : *(mp(il,i)-mp(il,i+1))*xt(ixt_poubelle,il,i) ! fxt(ixt_ddft,il,i)=fxt(ixt_ddft,il,i)+conversion(iiso) ! fxt(ixt_poubelle,il,i)=fxt(ixt_poubelle,il,i) ! : -conversion(iiso) ! Pb: quand on discretise, dqp/dt n'est pas verifee numeriquement. ! on se retrouve donc avec des d Ye/dt differents de 0 meme si ye=0 ( on ! note X les molecules poubelles et Y les molecules ddfts). ! Solution alternative: Dans le cas entrainant, Ye ne varie que par ! ascendance compensatoire des ddfts et par perte de Ye vers le ddft. On ! calcule donc ce terme directement avec schema amont: ! ajout deja de l'evap DO ixt = 1+niso,ntraciso fxt(ixt,il,i)=fxt(ixt,il,i) & +0.5*sigd*(xtevap(ixt,il,i)+xtevap(ixt,il,i+1)) enddo !do ixt = 1+niso,ntraciso ! ajout du terme des ddfts sensi stricto ! WRITE(*,*) 'tmp cv3_yield 4165: i,il=',i,il IF (option_traceurs.EQ.6) THEN DO iiso = 1, niso ixt_ddft=itZonIso(izone_ddft,iiso) IF (mp(il,i).gt.mp(il,i+1)) THEN fxtYe(iiso)=0.01*grav*dpinv*mp(il,i) & *(xt(ixt_ddft,il,i-1)-xt(ixt_ddft,il,i)) else !if (mp(il,i).gt.mp(il,i+1)) THEN fxtYe(iiso)=0.01*grav*dpinv*(mp(il,i) & *xt(ixt_ddft,il,i-1)-mp(il,i+1)*xt(ixt_ddft,il,i) & +(mp(il,i+1)-mp(il,i))*xtp(ixt_ddft,il,i)) endif !if (mp(il,i).gt.mp(il,i+1)) THEN fxtqe(iiso)=0.01*grav*dpinv* & (mp(il,i+1)*(xtp(iiso,il,i+1)-xt(iiso,il,i)) & -mp(il,i)*(xtp(iiso,il,i)-xt(iiso,il,i-1))) ixt_poubelle=itZonIso(izone_poubelle,iiso) fxt(ixt_ddft,il,i)=fxt(ixt_ddft,il,i)+fxtYe(iiso) fxt(ixt_poubelle,il,i)=fxt(ixt_poubelle,il,i) & +fxtqe(iiso)-fxtYe(iiso) enddo !do iiso = 1, niso else !if (option_traceurs.EQ.6) THEN IF (mp(il,i).gt.mp(il,i+1)) THEN ! cas entrainant: faire attention DO iiso = 1, niso fxtqe(iiso)=0.01*grav*dpinv* & (mp(il,i+1)*(xtp(iiso,il,i+1)-xt(iiso,il,i)) & -mp(il,i)*(xtp(iiso,il,i)-xt(iiso,il,i-1))) ixt_ddft=itZonIso(izone_ddft,iiso) fxtYe(iiso)=0.01*grav*dpinv*mp(il,i) & *(xt(ixt_ddft,il,i-1)-xt(ixt_ddft,il,i)) fxt(ixt_ddft,il,i)=fxt(ixt_ddft,il,i)+fxtYe(iiso) ixt_revap=itZonIso(izone_revap,iiso) fxt_revap(iiso)=0.01*grav*dpinv*(mp(il,i+1)* & (xtp(ixt_revap,il,i+1)-xt(ixt_revap,il,i)) & -mp(il,i)*(xtp(ixt_revap,il,i)-xt(ixt_revap,il,i-1))) fxt(ixt_revap,il,i)=fxt(ixt_revap,il,i) & +fxt_revap(iiso) fxtXe(iiso)=fxtqe(iiso)-fxtYe(iiso)-fxt_revap(iiso) Xe(iiso)=xt(iiso,il,i) & -xt(ixt_ddft,il,i)-xt(ixt_revap,il,i) IF (Xe(iiso).gt.ridicule) THEN DO izone=1,nzone IF ((izone.NE.izone_revap).AND. & (izone.NE.izone_ddft)) THEN ixt=itZonIso(izone,iiso) fxt(ixt,il,i)=fxt(ixt,il,i) & +xt(ixt,il,i)/Xe(iiso)*fxtXe(iiso) endif !if ((izone.NE.izone_revap).AND. enddo !do izone=1,nzone #ifdef ISOVERIF ! WRITE(*,*) 'iiso=',iiso ! WRITE(*,*) 'fxtqe=',fxtqe(iiso) ! WRITE(*,*) 'fxtYe=',fxtYe(iiso) ! WRITE(*,*) 'fxt_revap=',fxt_revap(iiso) ! WRITE(*,*) 'fxtXe=',fxtXe(iiso) ! WRITE(*,*) 'Xe=',Xe(iiso) ! WRITE(*,*) 'xt=',xt(:,il,i) CALL iso_verif_traceur_justmass(fxt(1,il,i), & 'cv30_routine 4646') #endif else !if (abs(dXe).gt.ridicule) THEN ! dans ce cas, fxtXe doit etre faible #ifdef ISOVERIF IF (delt*fxtXe(iiso).gt.ridicule) THEN WRITE(*,*) 'cv30_routines 6563: delt*fxtXe(iiso)=', & delt*fxtXe(iiso) stop endif #endif DO izone=1,nzone IF ((izone.NE.izone_revap).AND. & (izone.NE.izone_ddft)) THEN ixt=itZonIso(izone,iiso) IF (izone.EQ.izone_poubelle) THEN fxt(ixt,il,i)=fxt(ixt,il,i)+fxtXe(iiso) else !if (izone.EQ.izone_poubelle) THEN ! pas de tendance pour ce tag la endif !if (izone.EQ.izone_poubelle) THEN endif !if ((izone.NE.izone_revap).AND. enddo !do izone=1,nzone #ifdef ISOVERIF CALL iso_verif_traceur_justmass(fxt(1,il,i), & 'cv30_routine 4671') #endif endif !if (abs(dXe).gt.ridicule) THEN enddo !do iiso = 1, niso else !if (mp(il,i).gt.mp(il,i+1)) THEN ! cas detrainant: pas de problemes DO ixt=1+niso,ntraciso fxt(ixt,il,i)=fxt(ixt,il,i) & +0.01*grav*(mp(il,i+1)*(xtp(ixt,il,i+1)-xt(ixt,il,i)) & -mp(il,i)*(xtp(ixt,il,i)-xt(ixt,il,i-1)))*dpinv enddo !do ixt=1+niso,ntraciso #ifdef ISOVERIF CALL iso_verif_traceur_justmass(fxt(1,il,i), & 'cv30_routine 4685') #endif endif !if (mp(il,i).gt.mp(il,i+1)) THEN endif !if (option_traceurs.EQ.6) THEN ! WRITE(*,*) 'delt*conversion=',delt*conversion(iso_eau) ! WRITE(*,*) 'delt*fxtYe=',delt*fxtYe(iso_eau) ! WRITE(*,*) 'delt*fxtqe=',delt*fxtqe(iso_eau) endif ! if ((option_traceurs.NE.6).AND.(option_traceurs.NE.19)) THEN #endif ! cam verif #ifdef ISOVERIF DO ixt=1,niso CALL iso_verif_noNAN(fxt(ixt,il,i),'cv30_routines 3496') enddo #endif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,i), & fr(il,i),'cv30_routines 3493',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN IF (1.EQ.0) THEN IF ((iso_HDO.gt.0).AND.(delt*fr(il,i).gt.ridicule)) THEN IF (iso_verif_aberrant_enc_nostop( & fxt(iso_HDO,il,i)/fr(il,i), & 'cv30_yield 3662').EQ.1) THEN WRITE(*,*) 'il,i,icb(il),inb(il)=',il,i,icb(il),inb(il) WRITE(*,*) 'fr(il,i),delt=',fr(il,i),delt #ifdef DIAGISO IF (fq_ddft(il,i).NE.0.0) THEN WRITE(*,*) 'fq_ddft,deltaD=',fq_ddft(il,i),deltaD( & fxt_ddft(iso_HDO,il,i)/fq_ddft(il,i)) endif !if (fq_ddft(il,i).NE.0.0) THEN IF (fq_evapprecip(il,i).NE.0.0) THEN WRITE(*,*) 'fq_evapprecip,deltaD=',fq_evapprecip(il,i), & deltaD(fxt_evapprecip(iso_HDO,il,i) & /fq_evapprecip(il,i)) endif !if (fq_evapprecip(il,i).NE.0.0) THEN #endif endif !if (iso_verif_aberrant_enc_nostop( endif !if (iso_HDO.gt.0) THEN endif !if (1.EQ.0) THEN IF ((iso_HDO.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN IF (iso_verif_aberrant_enc_nostop((xt(iso_HDO,il,i) & +delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)),'cv30_yield 3757, ddfts') & .EQ.1) THEN WRITE(*,*) 'i,il,q,deltaD=',i,il,rr(il,i),deltaD( & xt(iso_HDO,il,i)/rr(il,i)) WRITE(*,*) 'i,il,fr,deltaD=',i,il,fr(il,i),deltaD( & fxt(iso_HDO,il,i)/fr(il,i)) stop endif ! if (iso_verif_aberrant_enc_nostop endif !if (iso_HDO.gt.0) THEN IF ((iso_HDO.gt.0).AND.(iso_O18.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*fxt(iso_HDO,il,i))/(rr(il,i)+delt*fr(il,i)), & (xt(iso_O18,il,i)+delt*fxt(iso_O18,il,i))/(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 5250,O18, ddfts') endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC ! WRITE(*,*) 'tmp cv3_yield 4224: i,il=',i,il CALL iso_verif_traceur_justmass(fxt(1,il,i),'cv30_routine 4107') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,i)+delt*fxt(ixt,il,i) enddo IF (iso_verif_tracpos_choix_nostop(xtnew, & 'cv30_yield 4221',1e-5).EQ.1) THEN WRITE(*,*) 'delt*fxt(,il,i)=',delt*fxt(1:ntraciso:2,il,i) WRITE(*,*) 'delt*fxt(,il,i)=',delt*fxt(:,il,i) WRITE(*,*) 'xt(,il,i)=',xt(:,il,i) WRITE(*,*) 'delt,sigd,grav,dpinv=',delt,sigd,grav,dpinv WRITE(*,*) 'xtevap(,il,i)=',xtevap(:,il,i) WRITE(*,*) 'xtevap(,il,i+1)=',xtevap(:,il,i+1) WRITE(*,*) 'mp(il,i+1),mp(il,i)=',mp(il,i+1),mp(il,i) WRITE(*,*) 'xtp(,il,i)=',xtp(:,il,i) WRITE(*,*) 'xtp(,il,i+1)=',xtp(:,il,i+1) WRITE(*,*) 'xt(,il,i)=',xt(:,il,i) WRITE(*,*) 'xt(,il,i-1)=',xt(:,il,i-1) ! rappel: fxt(ixt,il,i)=fxt(ixt,il,i) ! 0.5*sigd*(xtevap(ixt,il,i)+xtevap(ixt,il,i+1)) ! : +0.01*grav*(mp(il,i+1)*(xtp(ixt,il,i+1)-xt(ixt,il,i)) ! : -mp(il,i)*(xtp(ixt,il,i)-xt(ixt,il,i-1)))*dpinv ! stop endif #endif #endif #endif ELSE ! cvflag_grav fr(il, i) = fr(il, i) + 0.5*sigd*(evap(il,i)+evap(il,i+1)) + & 0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)*(rp(il,i)-rr(il, & i-1)))*dpinv fu(il, i) = fu(il, i) + 0.1*(mp(il,i+1)*(up(il,i+1)-u(il, & i))-mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv fv(il, i) = fv(il, i) + 0.1*(mp(il,i+1)*(vp(il,i+1)-v(il, & i))-mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv #ifdef ISO DO ixt = 1, ntraciso fxt(ixt,il,i)=fxt(ixt,il,i) & +0.5*sigd*(xtevap(ixt,il,i)+xtevap(ixt,il,i+1)) & +0.1*(mp(il,i+1)*(xtp(ixt,il,i+1)-xt(ixt,il,i)) & -mp(il,i)*(xtp(ixt,il,i)-xt(ixt,il,i-1)))*dpinv enddo ! ixt=1,niso #ifdef ISOTRAC IF (option_traceurs.NE.6) THEN ! facile: on fait comme l'eau DO ixt = 1+niso,ntraciso fxt(ixt,il,i)=fxt(ixt,il,i) & +0.5*sigd*(xtevap(ixt,il,i)+xtevap(ixt,il,i+1)) & +0.01*grav*(mp(il,i+1)*(xtp(ixt,il,i+1)-xt(ixt,il,i)) & -mp(il,i)*(xtp(ixt,il,i)-xt(ixt,il,i-1)))*dpinv enddo !do ixt = 1+niso,ntraciso else !if (option_traceurs.NE.6) THEN ! taggage des ddfts: voir blabla + haut DO ixt = 1+niso,ntraciso fxt(ixt,il,i)=fxt(ixt,il,i) & +0.5*sigd*(xtevap(ixt,il,i)+xtevap(ixt,il,i+1)) enddo !do ixt = 1+niso,ntraciso ! WRITE(*,*) 'tmp cv3_yield 4165: i,il=',i,il ! ixt_poubelle=itZonIso(izone_poubelle,iso_eau) ! ixt_ddft=itZonIso(izone_ddft,iso_eau) ! WRITE(*,*) 'delt*fxt(ixt_poubelle,il,i)=', ! : delt*fxt(ixt_poubelle,il,i) ! WRITE(*,*) 'delt*fxt(ixt_ddft,il,i)=',delt*fxt(ixt_ddft,il,i) ! WRITE(*,*) 'xt(iso_eau,il,i)=',xt(iso_eau,il,i) DO iiso = 1, niso ixt_poubelle=itZonIso(izone_poubelle,iiso) ixt_ddft=itZonIso(izone_ddft,iiso) IF (mp(il,i).gt.mp(il,i+1)) THEN fxtYe(iiso)=0.01*grav*dpinv*mp(il,i) & *(xt(ixt_ddft,il,i-1)-xt(ixt_ddft,il,i)) else !if (mp(il,i).gt.mp(il,i+1)) THEN fxtYe(iiso)=0.01*grav*dpinv*(mp(il,i) & *xt(ixt_ddft,il,i-1)-mp(il,i+1)*xt(ixt_ddft,il,i) & +(mp(il,i+1)-mp(il,i))*xtp(ixt_ddft,il,i)) endif !if (mp(il,i).gt.mp(il,i+1)) THEN fxtqe(iiso)=0.01*grav*dpinv* & (mp(il,i+1)*(xtp(iiso,il,i+1)-xt(iiso,il,i)) & -mp(il,i)*(xtp(iiso,il,i)-xt(iiso,il,i-1))) fxt(ixt_ddft,il,i)=fxt(ixt_ddft,il,i)+fxtYe(iiso) fxt(ixt_poubelle,il,i)=fxt(ixt_poubelle,il,i) & +fxtqe(iiso)-fxtYe(iiso) enddo !do iiso = 1, niso ! WRITE(*,*) 'delt*conversion=',delt*conversion(iso_eau) ! WRITE(*,*) 'delt*fxtYe=',delt*fxtYe(iso_eau) ! WRITE(*,*) 'delt*fxtqe=',delt*fxtqe(iso_eau) endif !if (option_traceurs.EQ.6) THEN #endif #ifdef DIAGISO fq_evapprecip(il,i)=fq_evapprecip(il,i) & +0.5*sigd*(evap(il,i)+evap(il,i+1)) fq_ddft(il,i)=fq_ddft(il,i) & +0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i) & *(rp(il,i)-rr(il,i-1)))*dpinv DO ixt = 1, niso fxt_evapprecip(ixt,il,i)=fxt_evapprecip(ixt,il,i) & +0.5*sigd*(xtevap(ixt,il,i)+xtevap(ixt,il,i+1)) fxt_ddft(ixt,il,i)=fxt_ddft(ixt,il,i) & +0.1*(mp(il,i+1)*(xtp(ixt,il,i+1)-xt(ixt,il,i)) & -mp(il,i)*(xtp(ixt,il,i)-xt(ixt,il,i-1)))*dpinv enddo ! ixt=1,niso #endif ! cam verif #ifdef ISOVERIF DO ixt=1,niso CALL iso_verif_noNaN(fxt(ixt,il,i),'cv30_yield 5083') enddo #endif #ifdef ISOVERIF IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,i), & fr(il,i),'cv30_routines 3522',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN IF ((iso_HDO.gt.0).AND.(delt*fr(il,i).gt.ridicule)) THEN IF (iso_verif_aberrant_enc_nostop( & fxt(iso_HDO,il,i)/fr(il,i), & 'cv30_yield 3690').EQ.1) THEN WRITE(*,*) 'i,icb(il),inb(il)=',i,icb(il),inb(il) stop endif endif !if (iso_HDO.gt.0) THEN IF ((iso_HDO.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre((xt(iso_HDO,il,i) & +delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)),'cv30_yield 3757b, ddfts') endif !if (iso_HDO.gt.0) THEN IF ((iso_HDO.gt.0).AND.(iso_O18.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*fxt(iso_HDO,il,i))/(rr(il,i)+delt*fr(il,i)), & (xt(iso_O18,il,i)+delt*fxt(iso_O18,il,i))/(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 3757b,O18, ddfts') endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,i),'cv30_routine 4172') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,1)+delt*fxt(ixt,il,1) enddo IF (iso_verif_tracpos_choix_nostop(xtnew,'cv30_yield 4295',1e-5) & .EQ.1) THEN WRITE(*,*) 'il,i=',il,i endif ! CALL iso_verif_tracpos_choix(xtnew,'cv30_yield 4295',1e-5) #endif #endif ! end cam verif #endif END IF ! cvflag_grav END IF ! i END DO ! sb: interface with the cloud parameterization: ! cld DO k = i + 1, nl DO il = 1, ncum IF (k<=inb(il) .AND. i<=inb(il)) THEN ! cld ! (saturated downdrafts resulting from mixing) ! cld qcond(il, i) = qcond(il, i) + elij(il, k, i) ! cld nqcond(il, i) = nqcond(il, i) + 1. ! cld END IF ! cld END DO ! cld END DO ! cld ! (particular case: no detraining level is found) ! cld DO il = 1, ncum ! cld IF (i<=inb(il) .AND. nent(il,i)==0) THEN ! cld qcond(il, i) = qcond(il, i) + (1.-ep(il,i))*clw(il, i) ! cld nqcond(il, i) = nqcond(il, i) + 1. ! cld END IF ! cld END DO ! cld DO il = 1, ncum ! cld IF (i<=inb(il) .AND. nqcond(il,i)/=0.) THEN ! cld qcond(il, i) = qcond(il, i)/nqcond(il, i) ! cld END IF ! cld END DO ! do j=1,ntra ! do il=1,ncum ! if (i.le.inb(il)) THEN ! dpinv=1.0/(ph(il,i)-ph(il,i+1)) ! cpinv=1.0/cpn(il,i) ! if (cvflag_grav) THEN ! ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv ! : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) ! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) ! else ! ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv ! : *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j)) ! : -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j))) ! END IF ! END IF ! i ! enddo ! enddo 500 END DO ! *** move the detrainment at level inb down to level inb-1 *** ! *** in such a way as to preserve the vertically *** ! *** integrated enthalpy and water tendencies *** DO il = 1, ncum ! attention, on corrige un probleme C Risi IF (cvflag_grav) THEN ax = 0.01*grav*ment(il, inb(il), inb(il))*(hp(il,inb(il))-h(il,inb(il))+t(il, & inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il), & inb(il))))/(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1))) ft(il, inb(il)) = ft(il, inb(il)) - ax ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il,inb(il & ))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il, & inb(il)))) bx = 0.01*grav*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb( & il)))/(ph(il,inb(il))-ph(il,inb(il)+1)) fr(il, inb(il)) = fr(il, inb(il)) - bx fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+ & 1))/(ph(il,inb(il)-1)-ph(il,inb(il))) cx = 0.01*grav*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il & )))/(ph(il,inb(il))-ph(il,inb(il)+1)) fu(il, inb(il)) = fu(il, inb(il)) - cx fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+ & 1))/(ph(il,inb(il)-1)-ph(il,inb(il))) dx = 0.01*grav*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il & )))/(ph(il,inb(il))-ph(il,inb(il)+1)) fv(il, inb(il)) = fv(il, inb(il)) - dx fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+ & 1))/(ph(il,inb(il)-1)-ph(il,inb(il))) #ifdef ISO DO ixt = 1, ntraciso xtbx(ixt)=0.01*grav*ment(il,inb(il),inb(il)) & *(xtent(ixt,il,inb(il),inb(il)) & -xt(ixt,il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1)) fxt(ixt,il,inb(il))=fxt(ixt,il,inb(il))-xtbx(ixt) fxt(ixt,il,inb(il)-1)=fxt(ixt,il,inb(il)-1) & +xtbx(ixt)*(ph(il,inb(il))-ph(il,inb(il)+1)) & /(ph(il,inb(il)-1)-ph(il,inb(il))) enddo !do ixt = 1, niso #endif else !IF (cvflag_grav) ax = 0.1*ment(il, inb(il), inb(il))*(hp(il,inb(il))-h(il,inb(il))+t(il, & inb(il))*(cpv-cpd)*(rr(il,inb(il))-qent(il,inb(il), & inb(il))))/(cpn(il,inb(il))*(ph(il,inb(il))-ph(il,inb(il)+1))) ft(il, inb(il)) = ft(il, inb(il)) - ax ft(il, inb(il)-1) = ft(il, inb(il)-1) + ax*cpn(il, inb(il))*(ph(il,inb(il & ))-ph(il,inb(il)+1))/(cpn(il,inb(il)-1)*(ph(il,inb(il)-1)-ph(il, & inb(il)))) bx = 0.1*ment(il, inb(il), inb(il))*(qent(il,inb(il),inb(il))-rr(il,inb( & il)))/(ph(il,inb(il))-ph(il,inb(il)+1)) fr(il, inb(il)) = fr(il, inb(il)) - bx fr(il, inb(il)-1) = fr(il, inb(il)-1) + bx*(ph(il,inb(il))-ph(il,inb(il)+ & 1))/(ph(il,inb(il)-1)-ph(il,inb(il))) cx = 0.1*ment(il, inb(il), inb(il))*(uent(il,inb(il),inb(il))-u(il,inb(il & )))/(ph(il,inb(il))-ph(il,inb(il)+1)) fu(il, inb(il)) = fu(il, inb(il)) - cx fu(il, inb(il)-1) = fu(il, inb(il)-1) + cx*(ph(il,inb(il))-ph(il,inb(il)+ & 1))/(ph(il,inb(il)-1)-ph(il,inb(il))) dx = 0.1*ment(il, inb(il), inb(il))*(vent(il,inb(il),inb(il))-v(il,inb(il & )))/(ph(il,inb(il))-ph(il,inb(il)+1)) fv(il, inb(il)) = fv(il, inb(il)) - dx fv(il, inb(il)-1) = fv(il, inb(il)-1) + dx*(ph(il,inb(il))-ph(il,inb(il)+ & 1))/(ph(il,inb(il)-1)-ph(il,inb(il))) #ifdef ISO DO ixt = 1, ntraciso xtbx(ixt)=0.1*ment(il,inb(il),inb(il)) & *(xtent(ixt,il,inb(il),inb(il)) & -xt(ixt,il,inb(il)))/(ph(il,inb(il))-ph(il,inb(il)+1)) fxt(ixt,il,inb(il))=fxt(ixt,il,inb(il))-xtbx(ixt) fxt(ixt,il,inb(il)-1)=fxt(ixt,il,inb(il)-1) & +xtbx(ixt)*(ph(il,inb(il))-ph(il,inb(il)+1)) & /(ph(il,inb(il)-1)-ph(il,inb(il))) enddo !do ixt = 1, niso #endif endif !IF (cvflag_grav) #ifdef ISO #ifdef DIAGISO fq_detrainement(il,inb(il))=fq_detrainement(il,inb(il))-bx fq_detrainement(il,inb(il)-1)=fq_detrainement(il,inb(il)-1) & +bx*(ph(il,inb(il))-ph(il,inb(il)+1)) & /(ph(il,inb(il)-1)-ph(il,inb(il))) DO ixt = 1, niso fxt_detrainement(ixt,il,inb(il))= & fxt_detrainement(ixt,il,inb(il))-xtbx(ixt) fxt_detrainement(ixt,il,inb(il)-1)= & fxt_detrainement(ixt,il,inb(il)-1) & +xtbx(ixt)*(ph(il,inb(il))-ph(il,inb(il)+1)) & /(ph(il,inb(il)-1)-ph(il,inb(il))) enddo #endif ! cam verif #ifdef ISOVERIF DO ixt=1,niso CALL iso_verif_noNaN(fxt(ixt,il,inb(il)),'cv30_yield 5083') enddo IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,inb(il)), & fr(il,inb(il)),'cv30_routines 3638',errmax,errmaxrel) CALL iso_verif_egalite_choix(fxt(iso_eau,il,inb(il)-1), & fr(il,inb(il)-1),'cv30_routines 3640',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN IF ((iso_HDO.gt.0).AND. & (rr(il,inb(il))+delt*fr(il,inb(il)).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre( & (xt(iso_HDO,il,inb(il))+delt*fxt(iso_HDO,il,inb(il))) & /(rr(il,inb(il))+delt*fr(il,inb(il))), & 'cv30_yield 3921, en inb') IF (iso_O18.gt.0) THEN IF (iso_verif_O18_aberrant_nostop( & (xt(iso_HDO,il,inb(il))+delt*fxt(iso_HDO,il,inb(il))) & /(rr(il,inb(il))+delt*fr(il,inb(il))), & (xt(iso_O18,il,inb(il))+delt*fxt(iso_O18,il,inb(il))) & /(rr(il,inb(il))+delt*fr(il,inb(il))), & 'cv30_yield 3921O18, en inb').EQ.1) THEN WRITE(*,*) 'il,inb(il)=',il,inb(il) k_tmp=0.1*ment(il,inb(il),inb(il))/(ph(il,inb(il))-ph(il,inb(il)+1)) WRITE(*,*) 'fr,frprec=',fr(il,inb(il)),fr(il,inb(il))+bx WRITE(*,*) 'M,dt,k_tmp*dt=',k_tmp,delt,k_tmp*delt WRITE(*,*) 'q,qe=',rr(il,inb(il)),qent(il,inb(il),inb(il)) WRITE(*,*) 'r=',k_tmp*delt*qent(il,inb(il),inb(il))/rr(il,inb(il)) WRITE(*,*) 'deltaDR,Re=',deltaD(xt(iso_HDO,il,inb(il))/rr(il,inb(il))), & deltaD(xtent(iso_HDO,il,inb(il),inb(il))/qent(il,inb(il),inb(il))) WRITE(*,*) 'deltaO18R,Re=',deltaO(xt(iso_O18,il,inb(il))/rr(il,inb(il))), & deltaO(xtent(iso_O18,il,inb(il),inb(il))/qent(il,inb(il),inb(il))) stop endif !if (iso_verif_O18_aberrant_nostop endif !if (iso_O18.gt.0) THEN endif !if (iso_HDO.gt.0) THEN IF ((iso_HDO.gt.0).AND. & (rr(il,inb(il)-1)+delt*fr(il,inb(il)-1).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre( & (xt(iso_HDO,il,inb(il)-1) & +delt*fxt(iso_HDO,il,inb(il)-1)) & /(rr(il,inb(il)-1)+delt*fr(il,inb(il)-1)), & 'cv30_yield 3921b, en inb-1') IF (iso_O18.gt.0) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,inb(il)-1)+delt*fxt(iso_HDO,il,inb(il)-1)) & /(rr(il,inb(il)-1)+delt*fr(il,inb(il)-1)), & (xt(iso_O18,il,inb(il)-1)+delt*fxt(iso_O18,il,inb(il)-1)) & /(rr(il,inb(il)-1)+delt*fr(il,inb(il)-1)), & 'cv30_yield 3921cO18, en inb-1') endif endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,inb(il)-1), & 'cv30_routine 4364') CALL iso_verif_traceur_justmass(fxt(1,il,inb(il)), & 'cv30_routine 4364b') DO ixt=1,ntraciso xtnew(ixt)=xt(ixt,il,inb(il))+delt*fxt(ixt,il,inb(il)) enddo IF (iso_verif_tracpos_choix_nostop(xtnew,'cv30_yield 4492',1e-5) & .EQ.1) THEN WRITE(*,*) 'il,i=',il,i endif ! CALL iso_verif_tracpos_choix(xtnew,'cv30_yield 4492',1e-5) #endif #endif ! end cam verif #endif END DO ! do j=1,ntra ! do il=1,ncum ! ex=0.1*ment(il,inb(il),inb(il)) ! : *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j)) ! : /(ph(il,inb(il))-ph(il,inb(il)+1)) ! ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex ! ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j) ! : +ex*(ph(il,inb(il))-ph(il,inb(il)+1)) ! : /(ph(il,inb(il)-1)-ph(il,inb(il))) ! enddo ! enddo ! *** homoginize tendencies below cloud base *** DO il = 1, ncum asum(il) = 0.0 bsum(il) = 0.0 csum(il) = 0.0 dsum(il) = 0.0 #ifdef ISO frsum(il)=0.0 DO ixt=1,ntraciso fxtsum(ixt,il)=0.0 bxtsum(ixt,il)=0.0 enddo #endif END DO DO i = 1, nl DO il = 1, ncum IF (i<=(icb(il)-1)) THEN asum(il) = asum(il) + ft(il, i)*(ph(il,i)-ph(il,i+1)) bsum(il) = bsum(il) + fr(il, i)*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il, & 1)))*(ph(il,i)-ph(il,i+1)) csum(il) = csum(il) + (lv(il,i)+(cl-cpd)*(t(il,i)-t(il, & 1)))*(ph(il,i)-ph(il,i+1)) dsum(il) = dsum(il) + t(il, i)*(ph(il,i)-ph(il,i+1))/th(il, i) #ifdef ISO frsum(il)=frsum(il)+fr(il,i) DO ixt=1,ntraciso fxtsum(ixt,il)=fxtsum(ixt,il)+fxt(ixt,il,i) bxtsum(ixt,il)=bxtsum(ixt,il)+fxt(ixt,il,i) & *(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1))) & *(ph(il,i)-ph(il,i+1)) enddo #endif END IF END DO END DO ! !!! do 700 i=1,icb(il)-1 DO i = 1, nl DO il = 1, ncum IF (i<=(icb(il)-1)) THEN ft(il, i) = asum(il)*t(il, i)/(th(il,i)*dsum(il)) fr(il, i) = bsum(il)/csum(il) #ifdef ISO IF (abs(csum(il)).gt.0.0) THEN DO ixt=1,ntraciso fxt(ixt,il,i)=bxtsum(ixt,il)/csum(il) enddo else !if (frsum(il).gt.ridicule) THEN IF (abs(frsum(il)).gt.0.0) THEN DO ixt=1,ntraciso fxt(ixt,il,i)=fr(il,i)*fxtsum(ixt,il)/frsum(il) enddo else !if (abs(frsum(il)).gt.0.0) THEN IF (abs(fr(il,i))*delt.gt.ridicule) THEN WRITE(*,*) 'cv30_yield 4048: fr(il,i)=',fr(il,i) stop else !if (abs(fr(il,i))*delt.gt.ridicule) THEN DO ixt=1,ntraciso fxt(ixt,il,i)=0.0 enddo IF (iso_eau.gt.0) THEN fxt(iso_eau,il,i)=1.0 endif endif !if (abs(fr(il,i))*delt.gt.ridicule) THEN endif !if (abs(frsum(il)).gt.0.0) THEN endif !if (frsum(il).gt.0) THEN #endif END IF END DO END DO #ifdef ISO #ifdef ISOVERIF DO i=1,nl DO il=1,ncum DO ixt=1,ntraciso CALL iso_verif_noNAN(fxt(ixt,il,i),'cv30_yield 3826') enddo enddo enddo #endif #ifdef ISOVERIF DO i=1,nl ! WRITE(*,*) 'cv30_routines temp 3967: i=',i DO il=1,ncum ! WRITE(*,*) 'cv30_routines 3969: il=',il ! WRITE(*,*) 'cv30_routines temp 3967: il,i,inb(il),ncum=', ! : il,i,inb(il),ncum ! WRITE(*,*) 'cv30_routines 3974' IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,i), & fr(il,i),'cv30_yield 3830',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN ! WRITE(*,*) 'cv30_routines 3979' IF ((iso_HDO.gt.0).AND. & (delt*fr(il,i).gt.ridicule)) THEN IF (iso_verif_aberrant_enc_nostop( & fxt(iso_HDO,il,i)/fr(il,i), & 'cv30_yield 3834').EQ.1) THEN IF (fr(il,i).gt.ridicule*1e5) THEN WRITE(*,*) 'il,i,icb(il)=',il,i,icb(il) WRITE(*,*) 'frsum(il)=',frsum(il) WRITE(*,*) 'fr(il,i)=',fr(il,i) WRITE(*,*) 'csum(il)=',csum(il) WRITE(*,*) & 'deltaD(bxtsum(iso_HDO,il)/csum(il))=', & deltaD(bxtsum(iso_HDO,il)/csum(il)) ! stop endif ! WRITE(*,*) 'cv30_routines 3986: temporaire' endif !if (iso_verif_aberrant_enc_nostop endif !if (iso_HDO.gt.0) THEN IF ((iso_HDO.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN IF (iso_verif_aberrant_enc_nostop( & (xt(iso_HDO,il,i)+delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)),'cv30_yield 3921c, dans la CL') & .EQ.1) THEN WRITE(*,*) 'il,i,icb(il)=',il,i,icb(il) WRITE(*,*) 'frsum(il)=',frsum(il) WRITE(*,*) 'fr(il,i)=',fr(il,i) stop endif endif !if (iso_HDO.gt.0) THEN IF ((iso_HDO.gt.0).AND.(iso_O18.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*fxt(iso_HDO,il,i))/(rr(il,i)+delt*fr(il,i)), & (xt(iso_O18,il,i)+delt*fxt(iso_O18,il,i))/(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 3921d, dans la CL') endif !if (iso_HDO.gt.0) THEN #ifdef ISOTRAC CALL iso_verif_traceur_justmass(fxt(1,il,i), & 'cv30_routine 4523') #endif ! WRITE(*,*) 'cv30_routines 3994' enddo !do il=1,ncum ! WRITE(*,*) 'cv30_routine 3990: fin des il pour i=',i enddo !do i=1,nl ! WRITE(*,*) 'cv30_routine 3990: fin des verifs sur homogen' #endif #ifdef ISOVERIF ! verif finale des tendances: DO i=1,nl DO il=1,ncum IF (iso_eau.gt.0) THEN CALL iso_verif_egalite_choix(fxt(iso_eau,il,i), & fr(il,i),'cv30_yield 3830',errmax,errmaxrel) endif !if (iso_eau.gt.0) THEN IF ((iso_HDO.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_aberrant_encadre((xt(iso_HDO,il,i) & +delt*fxt(iso_HDO,il,i)) & /(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 5710a, final') endif !if (iso_HDO.gt.0) THEN IF ((iso_HDO.gt.0).AND.(iso_O18.gt.0).AND. & (rr(il,i)+delt*fr(il,i).gt.ridicule)) THEN CALL iso_verif_O18_aberrant( & (xt(iso_HDO,il,i)+delt*fxt(iso_HDO,il,i))/(rr(il,i)+delt*fr(il,i)), & (xt(iso_O18,il,i)+delt*fxt(iso_O18,il,i))/(rr(il,i)+delt*fr(il,i)), & 'cv30_yield 5710b, final') endif !if (iso_HDO.gt.0) THEN enddo !do il=1,ncum enddo !do i=1,nl #endif #endif ! *** reset counter and return *** DO il = 1, ncum sig(il, nd) = 2.0 END DO DO i = 1, nd DO il = 1, ncum upwd(il, i) = 0.0 dnwd(il, i) = 0.0 END DO END DO DO i = 1, nl DO il = 1, ncum dnwd0(il, i) = -mp(il, i) END DO END DO DO i = nl + 1, nd DO il = 1, ncum dnwd0(il, i) = 0. END DO END DO DO i = 1, nl DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il)) THEN upwd(il, i) = 0.0 dnwd(il, i) = 0.0 END IF END DO END DO DO i = 1, nl DO k = 1, nl DO il = 1, ncum up1(il, k, i) = 0.0 dn1(il, k, i) = 0.0 END DO END DO END DO DO i = 1, nl DO k = i, nl DO n = 1, i - 1 DO il = 1, ncum IF (i>=icb(il) .AND. i<=inb(il) .AND. k<=inb(il)) THEN up1(il, k, i) = up1(il, k, i) + ment(il, n, k) dn1(il, k, i) = dn1(il, k, i) - ment(il, k, n) END IF END DO END DO END DO END DO DO i = 2, nl DO k = i, nl DO il = 1, ncum ! test if (i.ge.icb(il).AND.i.le.inb(il).AND.k.le.inb(il)) ! THEN IF (i<=inb(il) .AND. k<=inb(il)) THEN upwd(il, i) = upwd(il, i) + m(il, k) + up1(il, k, i) dnwd(il, i) = dnwd(il, i) + dn1(il, k, i) END IF END DO END DO END DO ! !!! DO il=1,ncum ! !!! do i=icb(il),inb(il) ! !!! ! !!! upwd(il,i)=0.0 ! !!! dnwd(il,i)=0.0 ! !!! do k=i,inb(il) ! !!! up1=0.0 ! !!! dn1=0.0 ! !!! do n=1,i-1 ! !!! up1=up1+ment(il,n,k) ! !!! dn1=dn1-ment(il,k,n) ! !!! enddo ! !!! upwd(il,i)=upwd(il,i)+m(il,k)+up1 ! !!! dnwd(il,i)=dnwd(il,i)+dn1 ! !!! enddo ! !!! enddo ! !!! ! !!! ENDDO ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! determination de la variation de flux ascendant entre ! deux niveau non dilue mike ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc DO i = 1, nl DO il = 1, ncum mike(il, i) = m(il, i) END DO END DO DO i = nl + 1, nd DO il = 1, ncum mike(il, i) = 0. END DO END DO DO i = 1, nd DO il = 1, ncum ma(il, i) = 0 END DO END DO DO i = 1, nl DO j = i, nl DO il = 1, ncum ma(il, i) = ma(il, i) + m(il, j) END DO END DO END DO DO i = nl + 1, nd DO il = 1, ncum ma(il, i) = 0. END DO END DO DO i = 1, nl DO il = 1, ncum IF (i<=(icb(il)-1)) THEN ma(il, i) = 0 END IF END DO END DO ! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! icb represente de niveau ou se trouve la ! base du nuage , et inb le top du nuage ! ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc DO i = 1, nd DO il = 1, ncum mke(il, i) = upwd(il, i) + dnwd(il, i) END DO END DO DO i = 1, nd DO il = 1, ncum rdcp = (rrd*(1.-rr(il,i))-rr(il,i)*rrv)/(cpd*(1.-rr(il, & i))+rr(il,i)*cpv) tls(il, i) = t(il, i)*(1000.0/p(il,i))**rdcp tps(il, i) = tp(il, i) END DO END DO ! *** diagnose the in-cloud mixing ratio *** ! cld ! *** of condensed water *** ! cld ! cld DO i = 1, nd ! cld DO il = 1, ncum ! cld mac(il, i) = 0.0 ! cld wa(il, i) = 0.0 ! cld siga(il, i) = 0.0 ! cld sax(il, i) = 0.0 ! cld END DO ! cld END DO ! cld DO i = minorig, nl ! cld DO k = i + 1, nl + 1 ! cld DO il = 1, ncum ! cld IF (i<=inb(il) .AND. k<=(inb(il)+1)) THEN ! cld mac(il, i) = mac(il, i) + m(il, k) ! cld END IF ! cld END DO ! cld END DO ! cld END DO ! cld DO i = 1, nl ! cld DO j = 1, i ! cld DO il = 1, ncum ! cld IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld .AND. j>=icb(il)) THEN ! cld sax(il, i) = sax(il, i) + rrd*(tvp(il,j)-tv(il,j)) & ! cld *(ph(il,j)-ph(il,j+1))/p(il, j) ! cld END IF ! cld END DO ! cld END DO ! cld END DO ! cld DO i = 1, nl ! cld DO il = 1, ncum ! cld IF (i>=icb(il) .AND. i<=(inb(il)-1) & ! cld .AND. sax(il,i)>0.0) THEN ! cld wa(il, i) = sqrt(2.*sax(il,i)) ! cld END IF ! cld END DO ! cld END DO ! cld DO i = 1, nl ! cld DO il = 1, ncum ! cld IF (wa(il,i)>0.0) & ! cld siga(il, i) = mac(il, i)/wa(il, i) & ! cld *rrd*tvp(il, i)/p(il, i)/100./delta ! cld siga(il, i) = min(siga(il,i), 1.0) ! cld ! IM cf. FH IF (iflag_clw==0) THEN qcondc(il, i) = siga(il, i)*clw(il, i)*(1.-ep(il,i)) & ! cld +(1.-siga(il,i))*qcond(il, i) ! cld ELSE IF (iflag_clw==1) THEN qcondc(il, i) = qcond(il, i) ! cld END IF END DO ! cld END DO ! cld END SUBROUTINE cv30_yield !RomP >>> SUBROUTINE cv30_tracer(nloc, len, ncum, nd, na, ment, sij, da, phi, phi2, & d1a, dam, ep, vprecip, elij, clw, epmlmmm, eplamm, icb, inb) IMPLICIT NONE ! inputs: INTEGER ncum, nd, na, nloc, len REAL ment(nloc, na, na), sij(nloc, na, na) REAL clw(nloc, nd), elij(nloc, na, na) REAL ep(nloc, na) INTEGER icb(nloc), inb(nloc) REAL vprecip(nloc, nd+1) ! ouputs: REAL da(nloc, na), phi(nloc, na, na) REAL phi2(nloc, na, na) REAL d1a(nloc, na), dam(nloc, na) REAL epmlmmm(nloc, na, na), eplamm(nloc, na) ! variables pour tracer dans precip de l'AA et des mel ! local variables: INTEGER i, j, k, nam1 REAL epm(nloc, na, na) nam1=na-1 ! Introduced because ep is not defined for j=na ! variables d'Emanuel : du second indice au troisieme ! ---> tab(i,k,j) -> de l origine k a l arrivee j ! ment, sij, elij ! variables personnelles : du troisieme au second indice ! ---> tab(i,j,k) -> de k a j ! phi, phi2 ! initialisations DO j = 1, na DO i = 1, ncum da(i, j) = 0. d1a(i, j) = 0. dam(i, j) = 0. eplamm(i, j) = 0. END DO END DO DO k = 1, na DO j = 1, na DO i = 1, ncum epm(i, j, k) = 0. epmlmmm(i, j, k) = 0. phi(i, j, k) = 0. phi2(i, j, k) = 0. END DO END DO END DO ! fraction deau condensee dans les melanges convertie en precip : epm ! et eau condensee precipitee dans masse d'air sature : l_m*dM_m/dzdz.dzdz DO j = 1, nam1 DO k = 1, j - 1 DO i = 1, ncum IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN !jyg epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j) epm(i, j, k) = 1. - (1.-ep(i,j))*clw(i, j)/max(elij(i,k,j), 1.E-16) epm(i, j, k) = max(epm(i,j,k), 0.0) END IF END DO END DO END DO DO j = 1, nam1 DO k = 1, nam1 DO i = 1, ncum IF (k>=icb(i) .AND. k<=inb(i)) THEN eplamm(i, j) = eplamm(i, j) + ep(i, j)*clw(i, j)*ment(i, j, k)*(1.- & sij(i,j,k)) END IF END DO END DO END DO DO j = 1, nam1 DO k = 1, j - 1 DO i = 1, ncum IF (k>=icb(i) .AND. k<=inb(i) .AND. j<=inb(i)) THEN epmlmmm(i, j, k) = epm(i, j, k)*elij(i, k, j)*ment(i, k, j) END IF END DO END DO END DO ! matrices pour calculer la tendance des concentrations dans cvltr.F90 DO j = 1, nam1 DO k = 1, nam1 DO i = 1, ncum da(i, j) = da(i, j) + (1.-sij(i,k,j))*ment(i, k, j) phi(i, j, k) = sij(i, k, j)*ment(i, k, j) d1a(i, j) = d1a(i, j) + ment(i, k, j)*ep(i, k)*(1.-sij(i,k,j)) END DO END DO END DO DO j = 1, nam1 DO k = 1, j - 1 DO i = 1, ncum dam(i, j) = dam(i, j) + ment(i, k, j)*epm(i, j, k)*(1.-ep(i,k))*(1.- & sij(i,k,j)) phi2(i, j, k) = phi(i, j, k)*epm(i, j, k) END DO END DO END DO END SUBROUTINE cv30_tracer ! RomP <<< SUBROUTINE cv30_uncompress(nloc, len, ncum, nd, ntra, idcum, iflag, precip, & vprecip, evap, ep, sig, w0, ft, fq, fu, fv, ftra, inb, ma, upwd, dnwd, & dnwd0, qcondc, wd, cape, da, phi, mp, phi2, d1a, dam, sij, elij, clw, & epmlmmm, eplamm, wdtraina, wdtrainm,epmax_diag, iflag1, precip1, vprecip1, evap1, & ep1, sig1, w01, ft1, fq1, fu1, fv1, ftra1, inb1, ma1, upwd1, dnwd1, & dnwd01, qcondc1, wd1, cape1, da1, phi1, mp1, phi21, d1a1, dam1, sij1, & elij1, clw1, epmlmmm1, eplamm1, wdtraina1, wdtrainm1,epmax_diag1 & ! epmax_cape #ifdef ISO ,xtprecip,fxt,xtVPrecip,xtevap,xtclw,xtwdtraina & ,xtprecip1,fxt1,xtVPrecip1,xtevap1,xtclw1,xtwdtraina1 & #ifdef DIAGISO , water,xtwater,qp,xtp & , fq_detrainement,fq_ddft,fq_fluxmasse,fq_evapprecip & , fxt_detrainement,fxt_ddft,fxt_fluxmasse, fxt_evapprecip & , f_detrainement,q_detrainement,xt_detrainement & , water1,xtwater1,qp1,xtp1 & , fq_detrainement1,fq_ddft1,fq_fluxmasse1,fq_evapprecip1 & , fxt_detrainement1,fxt_ddft1,fxt_fluxmasse1, fxt_evapprecip1 & , f_detrainement1,q_detrainement1,xt_detrainement1 & #endif #endif ) #ifdef ISO USE infotrac_phy, ONLY: ntraciso=>ntiso #ifdef ISOVERIF USE isotopes_verif_mod, ONLY: Tmin_verif,iso_verif_aberrant, & iso_verif_egalite,iso_verif_egalite_choix_nostop,iso_verif_positif_nostop, & iso_verif_egalite_nostop,iso_verif_aberrant_nostop,deltaD,iso_verif_noNaN_nostop, & iso_verif_positif,iso_verif_egalite_vect2D #endif #endif IMPLICIT NONE ! inputs: INTEGER len, ncum, nd, ntra, nloc INTEGER idcum(nloc) INTEGER iflag(nloc) INTEGER inb(nloc) REAL precip(nloc) REAL vprecip(nloc, nd+1), evap(nloc, nd) REAL ep(nloc, nd) REAL sig(nloc, nd), w0(nloc, nd) REAL ft(nloc, nd), fq(nloc, nd), fu(nloc, nd), fv(nloc, nd) REAL ftra(nloc, nd, ntra) REAL ma(nloc, nd) REAL upwd(nloc, nd), dnwd(nloc, nd), dnwd0(nloc, nd) REAL qcondc(nloc, nd) REAL wd(nloc), cape(nloc) REAL da(nloc, nd), phi(nloc, nd, nd), mp(nloc, nd) REAL epmax_diag(nloc) ! epmax_cape ! RomP >>> REAL phi2(nloc, nd, nd) REAL d1a(nloc, nd), dam(nloc, nd) REAL wdtraina(nloc, nd), wdtrainm(nloc, nd) REAL sij(nloc, nd, nd) REAL elij(nloc, nd, nd), clw(nloc, nd) REAL epmlmmm(nloc, nd, nd), eplamm(nloc, nd) ! RomP <<< #ifdef ISO REAL xtprecip(ntraciso,nloc) REAL xtvprecip(ntraciso,nloc, nd+1), xtevap(ntraciso,nloc, nd) REAL fxt(ntraciso,nloc,nd) REAL xtclw(ntraciso,nloc,nd) REAL xtwdtraina(ntraciso,nloc, nd) #endif ! outputs: INTEGER iflag1(len) INTEGER inb1(len) REAL precip1(len) REAL vprecip1(len, nd+1), evap1(len, nd) !<<< RomP REAL ep1(len, nd) !<<< RomP REAL sig1(len, nd), w01(len, nd) REAL ft1(len, nd), fq1(len, nd), fu1(len, nd), fv1(len, nd) REAL ftra1(len, nd, ntra) REAL ma1(len, nd) REAL upwd1(len, nd), dnwd1(len, nd), dnwd01(len, nd) REAL qcondc1(nloc, nd) REAL wd1(nloc), cape1(nloc) REAL da1(nloc, nd), phi1(nloc, nd, nd), mp1(nloc, nd) REAL epmax_diag1(len) ! epmax_cape ! RomP >>> REAL phi21(len, nd, nd) REAL d1a1(len, nd), dam1(len, nd) REAL wdtraina1(len, nd), wdtrainm1(len, nd) REAL sij1(len, nd, nd) REAL elij1(len, nd, nd), clw1(len, nd) REAL epmlmmm1(len, nd, nd), eplamm1(len, nd) ! RomP <<< #ifdef ISO REAL xtprecip1(ntraciso,len) REAL fxt1(ntraciso,len,nd) REAL xtVPrecip1(ntraciso,len,nd+1),xtevap1(ntraciso,len, nd) REAL xtwdtraina1(ntraciso,len, nd) REAL xtclw1(ntraciso,len, nd) #endif ! local variables: INTEGER i, k, j #ifdef ISO INTEGER ixt #endif #ifdef DIAGISO REAL water(nloc,nd) REAL xtwater(ntraciso,nloc,nd) REAL qp(nloc,nd),xtp(ntraciso,nloc,nd) REAL fq_detrainement(nloc,nd) REAL f_detrainement(nloc,nd) REAL q_detrainement(nloc,nd) REAL fq_ddft(nloc,nd) REAL fq_fluxmasse(nloc,nd) REAL fq_evapprecip(nloc,nd) REAL fxt_detrainement(ntraciso,nloc,nd) REAL xt_detrainement(ntraciso,nloc,nd) REAL fxt_ddft(ntraciso,nloc,nd) REAL fxt_fluxmasse(ntraciso,nloc,nd) REAL fxt_evapprecip(ntraciso,nloc,nd) REAL water1(len,nd) REAL xtwater1(ntraciso,len,nd) REAL qp1(len,nd),xtp1(ntraciso,len,nd) REAL fq_detrainement1(len,nd) REAL f_detrainement1(len,nd) REAL q_detrainement1(len,nd) REAL fq_ddft1(len,nd) REAL fq_fluxmasse1(len,nd) REAL fq_evapprecip1(len,nd) REAL fxt_detrainement1(ntraciso,len,nd) REAL xt_detrainement1(ntraciso,len,nd) REAL fxt_ddft1(ntraciso,len,nd) REAL fxt_fluxmasse1(ntraciso,len,nd) REAL fxt_evapprecip1(ntraciso,len,nd) #endif #ifdef ISOVERIF WRITE(*,*) 'cv30_routines 4293: entree dans cv3_uncompress' #endif DO i = 1, ncum precip1(idcum(i)) = precip(i) iflag1(idcum(i)) = iflag(i) wd1(idcum(i)) = wd(i) inb1(idcum(i)) = inb(i) cape1(idcum(i)) = cape(i) epmax_diag1(idcum(i))=epmax_diag(i) ! epmax_cape #ifdef ISO DO ixt = 1, ntraciso xtprecip1(ixt,idcum(i))=xtprecip(ixt,i) enddo #endif END DO DO k = 1, nl DO i = 1, ncum vprecip1(idcum(i), k) = vprecip(i, k) evap1(idcum(i), k) = evap(i, k) !<<< RomP sig1(idcum(i), k) = sig(i, k) w01(idcum(i), k) = w0(i, k) ft1(idcum(i), k) = ft(i, k) fq1(idcum(i), k) = fq(i, k) fu1(idcum(i), k) = fu(i, k) fv1(idcum(i), k) = fv(i, k) ma1(idcum(i), k) = ma(i, k) upwd1(idcum(i), k) = upwd(i, k) dnwd1(idcum(i), k) = dnwd(i, k) dnwd01(idcum(i), k) = dnwd0(i, k) qcondc1(idcum(i), k) = qcondc(i, k) da1(idcum(i), k) = da(i, k) mp1(idcum(i), k) = mp(i, k) ! RomP >>> ep1(idcum(i), k) = ep(i, k) d1a1(idcum(i), k) = d1a(i, k) dam1(idcum(i), k) = dam(i, k) clw1(idcum(i), k) = clw(i, k) eplamm1(idcum(i), k) = eplamm(i, k) wdtraina1(idcum(i), k) = wdtraina(i, k) wdtrainm1(idcum(i), k) = wdtrainm(i, k) ! RomP <<< #ifdef ISO DO ixt = 1, ntraciso fxt1(ixt,idcum(i),k)=fxt(ixt,i,k) xtVPrecip1(ixt,idcum(i),k)=xtVPrecip(ixt,i,k) xtevap1(ixt,idcum(i),k)=xtevap(ixt,i,k) xtwdtraina1(ixt,idcum(i),k)=xtwdtraina(ixt,i,k) xtclw1(ixt,idcum(i),k)=xtclw(ixt,i,k) enddo #endif END DO END DO DO i = 1, ncum sig1(idcum(i), nd) = sig(i, nd) END DO #ifdef ISO #ifdef DIAGISO DO k=1,nl DO i=1,ncum water1(idcum(i),k)=water(i,k) qp1(idcum(i),k)=qp(i,k) evap1(idcum(i),k)=evap(i,k) fq_detrainement1(idcum(i),k)=fq_detrainement(i,k) f_detrainement1(idcum(i),k)=f_detrainement(i,k) q_detrainement1(idcum(i),k)=q_detrainement(i,k) fq_ddft1(idcum(i),k)=fq_ddft(i,k) fq_evapprecip1(idcum(i),k)=fq_evapprecip(i,k) fq_fluxmasse1(idcum(i),k)=fq_fluxmasse(i,k) DO ixt = 1, ntraciso xtwater1(ixt,idcum(i),k)=xtwater(ixt,i,k) xtp1(ixt,idcum(i),k)=xtp(ixt,i,k) fxt_detrainement1(ixt,idcum(i),k)=fxt_detrainement(ixt,i,k) xt_detrainement1(ixt,idcum(i),k)=xt_detrainement(ixt,i,k) fxt_ddft1(ixt,idcum(i),k)=fxt_ddft(ixt,i,k) fxt_fluxmasse1(ixt,idcum(i),k)=fxt_fluxmasse(ixt,i,k) fxt_evapprecip1(ixt,idcum(i),k)=fxt_evapprecip(ixt,i,k) enddo enddo enddo DO i=1,ncum epmax_diag1(idcum(i))=epmax_diag(i) enddo #endif #endif ! do 2100 j=1,ntra ! do 2110 k=1,nd ! oct3 ! do 2120 i=1,ncum ! ftra1(idcum(i),k,j)=ftra(i,k,j) ! 2120 continue ! 2110 continue ! 2100 continue DO j = 1, nd DO k = 1, nd DO i = 1, ncum sij1(idcum(i), k, j) = sij(i, k, j) phi1(idcum(i), k, j) = phi(i, k, j) phi21(idcum(i), k, j) = phi2(i, k, j) elij1(idcum(i), k, j) = elij(i, k, j) epmlmmm1(idcum(i), k, j) = epmlmmm(i, k, j) END DO END DO END DO END SUBROUTINE cv30_uncompress SUBROUTINE cv30_epmax_fn_cape(nloc,ncum,nd & ,cape,ep,hp,icb,inb,clw,nk,t,h,lv & ,epmax_diag) USE conema3_mod_h USE cvthermo_mod_h IMPLICIT NONE ! On fait varier epmax en fn de la cape ! Il faut donc recalculer ep, et hp qui a deja ete calcule et ! qui en depend ! Toutes les autres variables fn de ep sont calculees plus bas. ! inputs: INTEGER ncum, nd, nloc INTEGER icb(nloc), inb(nloc) REAL cape(nloc) REAL clw(nloc,nd),lv(nloc,nd),t(nloc,nd),h(nloc,nd) INTEGER nk(nloc) ! inouts: REAL ep(nloc,nd) REAL hp(nloc,nd) ! outputs ou local REAL epmax_diag(nloc) ! locals INTEGER i,k REAL hp_bak(nloc,nd) CHARACTER (LEN=20) :: modname='cv30_epmax_fn_cape' CHARACTER (LEN=80) :: abort_message ! on recalcule ep et hp IF (coef_epmax_cape.gt.1e-12) THEN DO i=1,ncum epmax_diag(i)=epmax-coef_epmax_cape*sqrt(cape(i)) DO k=1,nl ep(i,k)=ep(i,k)/epmax*epmax_diag(i) ep(i,k)=amax1(ep(i,k),0.0) ep(i,k)=amin1(ep(i,k),epmax_diag(i)) enddo enddo ! On recalcule hp: DO k=1,nl DO i=1,ncum hp_bak(i,k)=hp(i,k) enddo enddo DO k=1,nlp DO i=1,ncum hp(i,k)=h(i,k) enddo enddo DO k=minorig+1,nl DO i=1,ncum IF((k.ge.icb(i)).AND.(k.le.inb(i)))THEN hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k) endif enddo enddo !do k=minorig+1,n ! WRITE(*,*) 'cv30_routines 6218: hp(1,20)=',hp(1,20) DO i=1,ncum DO k=1,nl IF (abs(hp_bak(i,k)-hp(i,k)).gt.0.01) THEN WRITE(*,*) 'i,k=',i,k WRITE(*,*) 'coef_epmax_cape=',coef_epmax_cape WRITE(*,*) 'epmax_diag(i)=',epmax_diag(i) WRITE(*,*) 'ep(i,k)=',ep(i,k) WRITE(*,*) 'hp(i,k)=',hp(i,k) WRITE(*,*) 'hp_bak(i,k)=',hp_bak(i,k) WRITE(*,*) 'h(i,k)=',h(i,k) WRITE(*,*) 'nk(i)=',nk(i) WRITE(*,*) 'h(i,nk(i))=',h(i,nk(i)) WRITE(*,*) 'lv(i,k)=',lv(i,k) WRITE(*,*) 't(i,k)=',t(i,k) WRITE(*,*) 'clw(i,k)=',clw(i,k) WRITE(*,*) 'cpd,cpv=',cpd,cpv CALL abort_physic(modname,abort_message,0) endif enddo !do k=1,nl enddo !do i=1,ncum endif !if (coef_epmax_cape.gt.1e-12) THEN END SUBROUTINE cv30_epmax_fn_cape END MODULE cv30_routines_mod