Index: trunk/LMDZ.MARS/libf/phymars/callkeys.h =================================================================== --- trunk/LMDZ.MARS/libf/phymars/callkeys.h (revision 148) +++ trunk/LMDZ.MARS/libf/phymars/callkeys.h (revision 161) @@ -12,5 +12,5 @@ & ,callg2d,linear,rayleigh,tracer,active,doubleq,submicron & & ,lifting,callddevil,scavenging,sedimentation,activice,water & - & ,caps,photochem + & ,caps,photochem,calltherm,outptherm COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche & @@ -23,5 +23,6 @@ & ,callstats,calleofdump & & ,callnirco2,callnlte,callthermos,callconduct, & - & calleuv,callmolvis,callmoldiff,thermochem,thermoswater + & calleuv,callmolvis,callmoldiff,thermochem,thermoswater & + & ,calltherm,outptherm Index: trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90 (revision 161) +++ trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90 (revision 161) @@ -0,0 +1,209 @@ +! +! AC 2011-01-05 +! + SUBROUTINE calltherm_interface (ngrid,nlayer,firstcall, & + & long,lati,zzlev,zzlay, & + & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, & + & pplay,pplev,pphi,nq,zpopsk, & + & pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,pbl_dtke) + + USE ioipsl_getincom + + implicit none +#include "callkeys.h" +!-------------------------------------------------------- +! Variables d'entree +!-------------------------------------------------------- + + INTEGER, INTENT(IN) :: ngrid,nlayer,nq + REAL, INTENT(IN) :: ptimestep + REAL, INTENT(IN) :: pplev(ngrid,nlayer+1),pplay(ngrid,nlayer) + REAL, INTENT(IN) :: pphi(ngrid,nlayer) + REAL, INTENT(IN) :: pu(ngrid,nlayer),pv(ngrid,nlayer) + REAL, INTENT(IN) :: pt(ngrid,nlayer),pq(ngrid,nlayer,nq) + REAL, INTENT(IN) :: zzlay(ngrid,nlayer) + REAL, INTENT(IN) :: zzlev(ngrid,nlayer+1) + LOGICAL, INTENT(IN) :: firstcall + REAL, INTENT(IN) :: pdu(ngrid,nlayer),pdv(ngrid,nlayer) + REAL, INTENT(IN) :: pdq(ngrid,nlayer,nq),pdt(ngrid,nlayer) + REAL, INTENT(IN) :: q2(ngrid,nlayer+1) + REAL, INTENT(IN) :: long(ngrid),lati(ngrid) + REAL, INTENT(IN) :: zpopsk(ngrid,nlayer) + +!-------------------------------------------------------- +! Variables de sortie (ou entree/sortie) +!-------------------------------------------------------- + + REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer) + REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq) + INTEGER lmax_th(ngrid) + REAL pbl_dtke(ngrid,nlayer+1) + +!-------------------------------------------------------- +! Variables du thermique +!-------------------------------------------------------- + REAL u_seri(ngrid,nlayer), v_seri(ngrid,nlayer) + REAL t_seri(ngrid,nlayer) + REAL d_t_ajs(ngrid,nlayer) + REAL d_u_ajs(ngrid,nlayer), d_q_ajs(ngrid,nlayer,nq) + REAL d_v_ajs(ngrid,nlayer) + REAL fm_therm(ngrid,nlayer+1), entr_therm(ngrid,nlayer) + REAL detr_therm(ngrid,nlayer) + REAL zw2(ngrid,nlayer+1) + REAL fraca(ngrid,nlayer+1) + REAL ztla(ngrid,nlayer) + REAL q_therm(ngrid,nlayer), pq_therm(ngrid,nlayer,nq) + REAL dq_therm(ngrid,nlayer), dq_thermdown(ngrid,nlayer) + REAL q2_therm(ngrid,nlayer), dq2_therm(ngrid,nlayer) + + LOGICAL qtransport_thermals,dtke_thermals + + INTEGER l,ig,iq + +! Variable de diagnostique : flux de chaleur vertical + + REAL heatFlux(ngrid,nlayer) + REAL heatFlux_down(ngrid,nlayer) + REAL buoyancyOut(ngrid,nlayer) + REAL buoyancyEst(ngrid,nlayer) + +!--------------------------------------------------------- +!--------------------------------------------------------- +! ********************************************************************** +! Thermique +! ********************************************************************** + +! Initialisation des sorties + + lmax_th(:)=1 + pdu_th(:,:)=0. + pdv_th(:,:)=0. + pdt_th(:,:)=0. + entr_therm(:,:)=0. + detr_therm(:,:)=0. + q2_therm(:,:)=0. + dq2_therm(:,:)=0. + dq_therm(:,:)=0. + dq_thermdown(:,:)=0. + ztla(:,:)=0. + pbl_dtke(:,:)=0. + fm_therm(:,:)=0. + zw2(:,:)=0. + fraca(:,:)=0. + if (tracer) then + pdq_th(:,:,:)=0. + end if + +! Dans le model terrestres, les seri sont des q+dq tendances déja cumulées. Il n'y a donc pas de +! cumulage à l'intérieur de la routine comme dans le model martien. On le fait ici : + + u_seri(:,:)=pu(:,:)+pdu(:,:)*ptimestep + v_seri(:,:)=pv(:,:)+pdv(:,:)*ptimestep + t_seri(:,:)=pt(:,:)+pdt(:,:)*ptimestep + + pq_therm(:,:,:)=0. + call getin("qtransport_thermals",qtransport_thermals) + if(qtransport_thermals) then + if(tracer) then + pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep + endif + endif + + d_t_ajs(:,:)=0. + d_u_ajs(:,:)=0. + d_v_ajs(:,:)=0. + d_q_ajs(:,:,:)=0. + heatFlux(:,:)=0. + heatFlux_down(:,:)=0. + buoyancyOut(:,:)=0. + buoyancyEst(:,:)=0. + + call getin("dtke_thermals",dtke_thermals) + if(dtke_thermals) then + + DO l=1,nlayer + q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1)) + ENDDO + endif + + CALL calltherm_mars(ngrid,nlayer,ptimestep,nq,zzlev,zzlay & + & ,pplay,pplev,pphi & + & ,u_seri,v_seri,t_seri,pq_therm, q2_therm & + & ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, dq2_therm & + & ,fm_therm,entr_therm,detr_therm & + & ,lmax_th & + & ,zw2,fraca & + & ,zpopsk,ztla,heatFlux,heatFlux_down & + & ,buoyancyOut,buoyancyEst) + + +! Accumulation des tendances. On n'accumule pas les quantités de traceurs car celle ci n'a pas du changer +! étant donné qu'on ne prends en compte que q_seri de la vap d'eau = 0 + +! INCREMENTATION : les d_u_ sont des tendances alors que les pdu sont des dérivees, attention ! + + pdu_th(:,:)=d_u_ajs(:,:)/ptimestep + pdv_th(:,:)=d_v_ajs(:,:)/ptimestep + pdt_th(:,:)=d_t_ajs(:,:)/ptimestep + if(qtransport_thermals) then + if(tracer) then + pdq_th(:,:,:)=d_q_ajs(:,:,:)/ptimestep + endif + endif + + + DO l=2,nlayer + pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))/ptimestep + ENDDO + + pbl_dtke(:,1)=0.5*dq2_therm(:,1)/ptimestep + pbl_dtke(:,nlayer+1)=0. +!! DIAGNOSTICS + + if(outptherm) then + if (ngrid .eq. 1) then + call WRITEDIAGFI(ngrid,'entr_therm','entrainement thermique',& + & 'kg/m-2',1,entr_therm) + call WRITEDIAGFI(ngrid,'detr_therm','detrainement thermique',& + & 'kg/m-2',1,detr_therm) + call WRITEDIAGFI(ngrid,'fm_therm','flux masse thermique',& + & 'kg/m-2',1,fm_therm) + call WRITEDIAGFI(ngrid,'zw2','vitesse verticale thermique',& + & 'm/s',1,zw2) + call WRITEDIAGFI(ngrid,'heatFlux_up','heatFlux_updraft',& + & 'SI',1,heatFlux) + call WRITEDIAGFI(ngrid,'heatFlux_down','heatFlux_downdraft',& + & 'SI',1,heatFlux_down) + call WRITEDIAGFI(ngrid,'fraca','fraction coverage',& + & 'percent',1,fraca) + call WRITEDIAGFI(ngrid,'buoyancyOut','buoyancyOut',& + & 'm.s-2',1,buoyancyOut) + call WRITEDIAGFI(ngrid,'buoyancyEst','buoyancyEst',& + & 'm.s-2',1,buoyancyEst) + call WRITEDIAGFI(ngrid,'d_t_th', & + & 'tendance temp TH','K',1,d_t_ajs) + call WRITEDIAGFI(ngrid,'d_q_ajs', & + & 'tendance q updraft','kg/kg',1,d_q_ajs(:,:,nq)) + call WRITEDIAGFI(ngrid,'d_q_tke', & + & 'tendance q updraft','q2/q2',1,pbl_dtke*ptimestep) + else + + call WRITEDIAGFI(ngrid,'entr_therm','entrainement thermique',& + & 'kg/m-2',3,entr_therm) + call WRITEDIAGFI(ngrid,'detr_therm','detrainement thermique',& + & 'kg/m-2',3,detr_therm) + call WRITEDIAGFI(ngrid,'fm_therm','flux masse thermique',& + & 'kg/m-2',3,fm_therm) + call WRITEDIAGFI(ngrid,'zw2','vitesse verticale thermique',& + & 'm/s',3,zw2) + call WRITEDIAGFI(ngrid,'heatFlux','heatFlux',& + & 'SI',3,heatFlux) + call WRITEDIAGFI(ngrid,'buoyancyOut','buoyancyOut',& + & 'SI',3,buoyancyOut) + call WRITEDIAGFI(ngrid,'d_t_th', & + & 'tendance temp TH','K',3,d_t_ajs) + + endif + endif + + END Index: trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90 (revision 161) +++ trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90 (revision 161) @@ -0,0 +1,199 @@ +! +! $Id: calltherm.F90 1428 2010-09-13 08:43:37Z fairhead $ +! + subroutine calltherm_mars(ngrid,nlayer,dtime,nq,zzlev,zzlay & + & ,pplay,paprs,pphi & + & ,u_seri,v_seri,t_seri,pq_therm,q2_therm & + & ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,dq2_therm & + & ,fm_therm,entr_therm,detr_therm,lmax,& + & zw2,fraca,zpopsk,ztla,heatFlux,heatFlux_down,& + & buoyancyOut,buoyancyEst) + + USE thermcell, only : nsplit_thermals,r_aspect_thermals + USE ioipsl_getincom + implicit none + + INTEGER, INTENT(IN) :: ngrid,nlayer + REAL dtime + LOGICAL logexpr0, logexpr2(ngrid,nlayer), logexpr1(ngrid) + REAL fact + INTEGER nbptspb,nq + + REAL, INTENT(IN) :: zzlay(ngrid,nlayer) + REAL, INTENT(IN) :: zzlev(ngrid,nlayer+1) + + REAL u_seri(ngrid,nlayer),v_seri(ngrid,nlayer) + REAL t_seri(ngrid,nlayer),pq_therm(ngrid,nlayer,nq) + REAL q2_therm(ngrid,nlayer) + REAL paprs(ngrid,nlayer+1) + REAL pplay(ngrid,nlayer) + REAL pphi(ngrid,nlayer) + real zlev(ngrid,nlayer+1) +!test: on sort lentr et a* pour alimenter KE + REAL zw2(ngrid,nlayer+1),fraca(ngrid,nlayer+1) + REAL zzw2(ngrid,nlayer+1) + +!FH Update Thermiques + REAL d_t_ajs(ngrid,nlayer), d_q_ajs(ngrid,nlayer,nq) + REAL d_u_ajs(ngrid,nlayer),d_v_ajs(ngrid,nlayer) + REAL dq2_therm(ngrid,nlayer), dq2_the(ngrid,nlayer) + real fm_therm(ngrid,nlayer+1) + real entr_therm(ngrid,nlayer),detr_therm(ngrid,nlayer) + +!******************************************************** +! declarations + LOGICAL flag_bidouille_stratocu + real fmc_therm(ngrid,nlayer+1) + real zqla(ngrid,nlayer) + real zqta(ngrid,nlayer) + real zpopsk(ngrid,nlayer) + real ztla(ngrid,nlayer) + real wmax_sec(ngrid) + real zmax_sec(ngrid) + real f_sec(ngrid) + real detrc_therm(ngrid,nlayer) + real zw_sec(ngrid,nlayer+1) + integer lmix_sec(ngrid) + integer lmax(ngrid) + +!nouvelles variables pour la convection +!RC + !on garde le zmax du pas de temps precedent +!******************************************************** + + +! variables locales + REAL d_t_the(ngrid,nlayer), d_q_the(ngrid,nlayer,nq) + REAL d_u_the(ngrid,nlayer),d_v_the(ngrid,nlayer) +! + integer isplit + real zfm_therm(ngrid,nlayer+1),zdt + real zentr_therm(ngrid,nlayer),zdetr_therm(ngrid,nlayer) + real heatFlux(ngrid,nlayer) + real heatFlux_down(ngrid,nlayer) + real buoyancyOut(ngrid,nlayer) + real buoyancyEst(ngrid,nlayer) + real zheatFlux(ngrid,nlayer) + real zheatFlux_down(ngrid,nlayer) + real zbuoyancyOut(ngrid,nlayer) + real zbuoyancyEst(ngrid,nlayer) + + character (len=20) :: modname='calltherm' + character (len=80) :: abort_message + + integer i,k + logical, save :: first=.true. + +! Modele du thermique +! =================== + + nsplit_thermals=20 + call getin("nsplit_thermals",nsplit_thermals) + + fm_therm(:,:)=0. + detr_therm(:,:)=0. + entr_therm(:,:)=0. + + heatFlux(:,:)=0. + heatFlux_down(:,:)=0. + buoyancyOut(:,:)=0. + buoyancyEst(:,:)=0. + + zw2(:,:)=0. + + zdt=dtime/REAL(nsplit_thermals) + + do isplit=1,nsplit_thermals + +! On reinitialise les flux de masse a zero pour le cumul en +! cas de splitting + + zfm_therm(:,:)=0. + zentr_therm(:,:)=0. + zdetr_therm(:,:)=0. + + zheatFlux(:,:)=0. + zheatFlux_down(:,:)=0. + zbuoyancyOut(:,:)=0. + zbuoyancyEst(:,:)=0. + + zzw2(:,:)=0. + + d_t_the(:,:)=0. + d_u_the(:,:)=0. + d_v_the(:,:)=0. + dq2_the(:,:)=0. + if (nq .ne. 0) then + d_q_the(:,:,:)=0. + endif + + CALL thermcell_main_mars(ngrid,nlayer,nq,zdt & + & ,pplay,paprs,pphi,zzlev,zzlay & + & ,u_seri,v_seri,t_seri,pq_therm,q2_therm & + & ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the & + & ,zfm_therm,zentr_therm,zdetr_therm,lmax & + & ,r_aspect_thermals & + & ,zzw2,fraca,zpopsk & + & ,ztla,zheatFlux,zheatFlux_down & + & ,zbuoyancyOut,zbuoyancyEst) + + fact=1./REAL(nsplit_thermals) +! transformation de la derivee en tendance + + d_t_the(:,:)=d_t_the(:,:)*dtime*fact + d_u_the(:,:)=d_u_the(:,:)*fact + d_v_the(:,:)=d_v_the(:,:)*fact + dq2_the(:,:)=dq2_the(:,:)*fact + + if (nq .ne. 0) then + d_q_the(:,:,:)=d_q_the(:,:,:)*fact + endif + + fm_therm(:,:)=fm_therm(:,:) & + & +zfm_therm(:,:)*fact + entr_therm(:,:)=entr_therm(:,:) & + & +zentr_therm(:,:)*fact + detr_therm(:,:)=detr_therm(:,:) & + & +zdetr_therm(:,:)*fact + + heatFlux(:,:)=heatFlux(:,:) & + & +zheatFlux(:,:)*fact + heatFlux_down(:,:)=heatFlux_down(:,:) & + +zheatFlux_down(:,:)*fact + buoyancyOut(:,:)=buoyancyOut(:,:) & + & +zbuoyancyOut(:,:)*fact + buoyancyEst(:,:)=buoyancyEst(:,:) & + & +zbuoyancyEst(:,:)*fact + + zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact + +! accumulation de la tendance + + d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:) + d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:) + d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:) + d_q_ajs(:,:,:)=d_q_ajs(:,:,:)+d_q_the(:,:,:) + dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:) +! incrementation des variables meteo + + t_seri(:,:) = t_seri(:,:) + d_t_the(:,:) + u_seri(:,:) = u_seri(:,:) + d_u_the(:,:) + v_seri(:,:) = v_seri(:,:) + d_v_the(:,:) + pq_therm(:,:,:) = pq_therm(:,:,:) + d_q_the(:,:,:) + q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:) + + enddo ! isplit + + +!**************************************************************** + +! do i=1,ngrid +! do k=1,nlayer +! if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0. +! print*,'youpi je sers a quelque chose !' +! enddo +! enddo + + return + + end Index: trunk/LMDZ.MARS/libf/phymars/convadj.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/convadj.F (revision 148) +++ trunk/LMDZ.MARS/libf/phymars/convadj.F (revision 161) @@ -1,4 +1,4 @@ subroutine convadj(ngrid,nlay,nq,ptimestep, - & pplay,pplev,ppopsk, + & pplay,pplev,ppopsk,lmax_th, & pu,pv,ph,pq, & pdufi,pdvfi,pdhfi,pdqfi, @@ -21,5 +21,6 @@ ! Modif. 2005 by F. Forget. ! Modif. 2010 by R. Wordsworth -! +! Modif. 2011 by A. Colaitis +! !================================================================== @@ -38,5 +39,5 @@ ! --------- - INTEGER,intent(in) :: ngrid,nlay + INTEGER,intent(in) :: ngrid,nlay,lmax_th(ngrid) REAL,intent(in) :: ptimestep REAL,intent(in) :: ph(ngrid,nlay) @@ -163,6 +164,7 @@ DO l=2,nlay DO ig=1,ngrid - IF(zhc(ig,l).LT.zhc(ig,l-1)) vtest(ig)=.true. - ENDDO + IF((zhc(ig,l).LT.zhc(ig,l-1)) $ + $ .AND. (l .GT. lmax_th(ig))) vtest(ig)=.true. + ENDDO ENDDO @@ -198,9 +200,9 @@ ! Test loop upwards on l2 - DO - l2 = l2 + 1 - IF (l2 .GT. nlay) EXIT - IF (zhc(i, l2) .LT. zhc(i, l2-1)) THEN - + DO + l2 = l2 + 1 + IF (l2 .GT. nlay) EXIT + IF ((zhc(i, l2) .LT. zhc(i, l2-1)).AND.(l2 .GT. lmax_th(i))) THEN + ! l2 is the highest level of the unstable column @@ -228,10 +230,10 @@ ! do we have to extend the column downwards? - down = .false. - IF (l1 .ne. 1) then !-- and then - IF (zhmc .lt. zhc(i, l1-1)) then - down = .true. - END IF - END IF + down = .false. + IF (l1 .ne. 1) then !-- and then + IF ((zhmc .lt. zhc(i, l1-1)).and.(l1.gt.lmax_th(i))) then + down = .true. + END IF + END IF ! this could be a problem... Index: trunk/LMDZ.MARS/libf/phymars/inifis.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/inifis.F (revision 148) +++ trunk/LMDZ.MARS/libf/phymars/inifis.F (revision 161) @@ -204,4 +204,14 @@ write(*,*) " calldifv = ",calldifv + write(*,*) "call thermals ?" + calltherm=.false. ! default value + call getin("calltherm",calltherm) + write(*,*) " calltherm = ",calltherm + + write(*,*) "output thermal diagnostics ?" + outptherm=.false. ! default value + call getin("outptherm",outptherm) + write(*,*) " outptherm = ",outptherm + write(*,*) "call convective adjustment ?" calladj=.true. ! default value @@ -209,4 +219,8 @@ write(*,*) " calladj = ",calladj + if (calltherm .and. (.not. calladj)) then + print*,'Convadj has to be activated when using thermals' + stop + endif write(*,*) "call CO2 condensation ?" Index: trunk/LMDZ.MARS/libf/phymars/physiq.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 148) +++ trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 161) @@ -6,5 +6,4 @@ $ pw, $ pdu,pdv,pdt,pdq,pdpsrf,tracerdyn) - IMPLICIT NONE @@ -301,4 +300,12 @@ REAL time_phys +c Variables from thermal + + REAL lmax_th_out(ngrid) + REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer) + REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq) + INTEGER lmax_th(ngrid) + REAL dtke_th(ngrid,nlayer+1) + REAL dummycol(ngrid) c======================================================================= @@ -515,4 +522,5 @@ c Radiative transfer c ------------------ + CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, @@ -599,7 +607,6 @@ c 4. Vertical diffusion (turbulent mixing): c ----------------------------------------- -c + IF (calldifv) THEN - DO ig=1,ngrid @@ -624,4 +631,5 @@ & zdqdif,zdqsdif) + DO l=1,nlayer DO ig=1,ngrid @@ -661,4 +669,43 @@ ENDIF ! of IF (calldifv) +c----------------------------------------------------------------------- +c TEST. Thermals : +c HIGHLY EXPERIMENTAL, BEWARE !! +c ----------------------------- + + if(calltherm) then + + call calltherm_interface(ngrid,nlayer,firstcall, + $ long,lati,zzlev,zzlay, + $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, + $ pplay,pplev,pphi,nq,zpopsk, + $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,dtke_th) + + DO l=1,nlayer + DO ig=1,ngrid + pdu(ig,l)=pdu(ig,l)+pdu_th(ig,l) + pdv(ig,l)=pdv(ig,l)+pdv_th(ig,l) + pdt(ig,l)=pdt(ig,l)+pdt_th(ig,l) + q2(ig,l)=q2(ig,l)+dtke_th(ig,l)*ptimestep + ENDDO + ENDDO + + DO ig=1,ngrid + q2(ig,nlayer+1)=q2(ig,nlayer+1)+dtke_th(ig,nlayer+1)*ptimestep + ENDDO + + if (tracer) then + DO iq=1,nq + DO l=1,nlayer + DO ig=1,ngrid + pdq(ig,l,iq)=pdq(ig,l,iq)+pdq_th(ig,l,iq) + ENDDO + ENDDO + ENDDO + endif + + else !of if calltherm + lmax_th(:)=0 + end if c----------------------------------------------------------------------- @@ -678,9 +725,10 @@ CALL convadj(ngrid,nlayer,nq,ptimestep, - $ pplay,pplev,zpopsk, + $ pplay,pplev,zpopsk,lmax_th, $ pu,pv,zh,pq, $ pdu,pdv,zdh,pdq, $ zduadj,zdvadj,zdhadj, $ zdqadj) + DO l=1,nlayer @@ -1144,5 +1192,4 @@ c which can later be used to make the statistic files of the run: c "stats") only possible in 3D runs ! - IF (callstats) THEN @@ -1261,4 +1308,6 @@ c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, c & emis) + call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay) + call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev) call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, & tsurf) @@ -1277,10 +1326,10 @@ & fluxtop_sw_tot) call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) -c call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) -c call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) -c call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw) + call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) + call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) + call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw) call WRITEDIAGFI(ngrid,"rho","density","none",3,rho) c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2) -c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh) + call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh) c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay) c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2, @@ -1290,4 +1339,8 @@ c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate', c & 'w.m-2',3,zdtlw) +c CALL WRITEDIAGFI(ngridmx,'tauTESap', +c & 'tau abs 825 cm-1', +c & '',2,tauTES) + c ---------------------------------------------------------- @@ -1402,4 +1455,27 @@ c ---------------------------------------------------------- +c Outputs of thermals +c ---------------------------------------------------------- + + +! call WRITEDIAGFI(ngrid,'dtke', +! & 'tendance tke thermiques','m**2/s**2', +! & 3,dtke_th) +! call WRITEDIAGFI(ngrid,'d_u_ajs', +! & 'tendance u thermiques','m/s', +! & 3,pdu_th*ptimestep) +! call WRITEDIAGFI(ngrid,'d_v_ajs', +! & 'tendance v thermiques','m/s', +! & 3,pdv_th*ptimestep) +! if (tracer) then +! if (nq .eq. 2) then +! call WRITEDIAGFI(ngrid,'deltaq_th', +! & 'delta q thermiques','kg/kg', +! & 3,ptimestep*pdq_th(:,:,2)) +! endif +! endif + + +c ---------------------------------------------------------- c Output in netcdf file "diagsoil.nc" for subterranean c variables (output every "ecritphy", as for writediagfi) @@ -1435,4 +1511,53 @@ c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1') +! THERMALS STUFF 1D + if(calltherm) then + + lmax_th_out(:)=real(lmax_th(:)) + + if (ngridmx .eq. 1) then + call WRITEDIAGFI(ngridmx,'lmax_th', + & 'hauteur du thermique','K', + & 0,lmax_th_out) + + else + call WRITEDIAGFI(ngridmx,'lmax_th', + & 'hauteur du thermique','K', + & 2,lmax_th_out) + + endif + + co2col(:)=0. + dummycol(:)=0. + if (tracer) then + do l=1,nlayermx + do ig=1,ngrid + co2col(ig)=co2col(ig)+ + & zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g + if (nqmx .gt. 1) then + dummycol(ig)=dummycol(ig)+ + & zq(ig,l,2)*(pplev(ig,l)-pplev(ig,l+1))/g + endif + enddo + enddo + + end if + call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' & + & ,'kg/m-2',0,co2col) + call WRITEDIAGFI(ngrid,'dummycol','integrated dummy mass' & + & ,'kg/m-2',0,dummycol) + endif + call WRITEDIAGFI(ngrid,'w','vertical velocity' & + & ,'m/s',1,pw) + call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2) + call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0, + & tsurf) + + call WRITEDIAGFI(ngrid,"fluxsurf_sw", + & "Solar radiative flux to surface","W.m-2",0, + & fluxsurf_sw_tot) + + call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay) + call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev) ! or output in diagfi.nc (for testphys1d) call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps) Index: trunk/LMDZ.MARS/libf/phymars/testphys1d.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/testphys1d.F (revision 148) +++ trunk/LMDZ.MARS/libf/phymars/testphys1d.F (revision 161) @@ -229,4 +229,22 @@ close(91) endif ! of if (txt.eq."co2") + ! Allow for an initial profile of argon + ! Can also be used to introduce a decaying tracer + ! in the 1D (TBD) to study thermals + if (txt.eq."ar") then + !look for a "profile_ar" input file + open(91,file='profile_ar',status='old', + & form='formatted',iostat=ierr) + if (ierr.eq.0) then + read(91,*) qsurf(iq) + do ilayer=1,nlayermx + read(91,*) q(ilayer,iq) + enddo + else + write(*,*) "No profile_ar file!" + endif + close(91) + endif ! of if (txt.eq."ar") + ! WATER VAPOUR if (txt.eq."h2o_vap") then Index: trunk/LMDZ.MARS/libf/phymars/thermcell.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/thermcell.F90 (revision 161) +++ trunk/LMDZ.MARS/libf/phymars/thermcell.F90 (revision 161) @@ -0,0 +1,108 @@ +MODULE thermcell +! USE ioipsl_getincom + IMPLICIT NONE + + INTEGER :: iflag_thermals + REAL :: r_aspect_thermals,l_mix_thermals,tau_thermals + INTEGER :: w2di_thermals + INTEGER :: nsplit_thermals + INTEGER :: isplit + INTEGER :: iflag_coupl,iflag_clos,iflag_wake + INTEGER :: iflag_thermals_ed,iflag_thermals_optflux + REAL :: RG,RD + REAL :: rlvtt,rcpd,rtt,r2es + REAL :: retv,r5les,r5ies,rkappa + REAL :: R4LES,R4IES,R3IES,R3LES + INTEGER :: klon,klev + INTEGER :: prt_level,lunout + real,allocatable :: rlatd(:) + real,allocatable :: rlond(:) + +! Thermodynamic constants. The [OK] means that the constant has been made +! compatible with the Martian model +! ---------------------------------- + +! RG : mars mean gravity field [OK] +! RD : dry air constant = 1000 R/Mair [OK] +! rlvtt : value of Lv(Tt) (vapo. latent heat at Tt) [dep. on air ?] +! rcpd : Cp of dry air = 7/2 RD for a perfect gas [OK,CHECK VALUE with gcm one] +!rcpd = r/(r/CP), avec r/CP = 1/4.4 pour LES, 0.256793 pour gcm (1./3.89419), et r=1000.R/Mair=RD +! rtt : triple point temperature Tt [OK] +! r2es : [Probably OK] +! retv = restt*RD/RV : [-] +! restt : saturation pressure at Tt es(Tt) = 611.14 Pa [dep. on air ?] +! RV : vapor constant = 1000 R/Mh2o [OK] +! r3les : [Probably OK] +! r3ies : [Probably OK] +! r4les : [Probably OK] +! r4ies : [Probably OK] +! r5ies = r3ies*(rtt-r4les) [Probably OK] +! rkappa = RD/RCPD = r/cp = 0.256793 [OK] +! FOEEW : vapeur d'eau saturante [OK] +! FOEDE : derive par rapport a la temperature [OK] + PARAMETER (RG=3.72,RD=191.182) + PARAMETER (rlvtt=2.5008E+06,rcpd=RD*3.89419) + PARAMETER (rtt=273.16,r2es=253.156) + PARAMETER (retv=1.41409,r5les=4097.93) + PARAMETER (r5ies=5807.81,rkappa=1./3.89419) + PARAMETER (R4IES=7.66,R4LES=35.86) + PARAMETER (R3IES=21.875,R3LES=17.269) + +! r_aspect_thermals : rapport d'aspect : parameter defined in FH's HDR +! length over height ratio of the thermals in the alimentation layer (check) +! it's value is advised to be set to 2 in HDR [OK] +! l_mix_thermals : mixing length => lambda : parameter for plume penetration in +! the atm (the plum size will decrease with sqrt(lambda) : USELESS +! w2di_thermals : ? +! tau_thermals : relaxation length of the thermals : the model implements the thermals +! tendancies with a typical relaxation time. Typical value is said to be 1800 s in the +! code however, it is found at 720 s in the .def files [OK] +! nsplit_thermals : subdivisions of the timestep in the calculations of the plume. +! Increase the calculation precision greatly, but also computation time +! typical value is found at 1 in .def (no split) [OK] +! prt_level : print level for the gcm printed outputs [OK] +! lunout : output channel for prt level [OK] +! iflag_thermals : choice of the thermals version. 15 and 16 are the newest ones +! and 16 activates "bidouille stratocu" which is required for now +! [OK] +! iflag_thermals_ed : choice of the plume version +! 8 is the normal version, 9 is the working version of AJ, 10 is CR +! also used for cases .ge.1 in thermcell_height [OK] +! iflag_coupl : USELESS ? [useless ?] +! iflag_clos : USELESS ? [useless ?] +! iflag_wake : USELESS ? [useless ?] + PARAMETER (r_aspect_thermals = 3.) + PARAMETER (l_mix_thermals = 30.) + PARAMETER (w2di_thermals = 1) +! PARAMETER (tau_thermals = 720.) +! PARAMETER (nsplit_thermals = 1) + PARAMETER (prt_level=0,lunout=6) + PARAMETER (iflag_thermals=20,iflag_thermals_ed=10) + PARAMETER (iflag_thermals_optflux=1,iflag_coupl=1) + PARAMETER (iflag_clos=2) + + CONTAINS + FUNCTION foede(PTARG,PDELARG,P5ARG,PQSARG,PCOARG) + IMPLICIT NONE + REAL :: foede + REAL, INTENT(IN) :: PTARG,PDELARG,P5ARG,PQSARG,PCOARG + foede = PQSARG*PCOARG*P5ARG & + & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG))**2 + END FUNCTION foede + + FUNCTION foeew(PTARG,PDELARG) + IMPLICIT NONE + REAL :: foeew + REAL, INTENT(IN) :: PTARG,PDELARG + foeew = EXP ( & + & (R3LES*(1.-PDELARG)+R3IES*PDELARG) * (PTARG-RTT) & + & / (PTARG-(R4LES*(1.-PDELARG)+R4IES*PDELARG)) ) + + END FUNCTION foeew + +! Variables which have been moved to .def and called independantly : +! nsplit_thermals, tau + +END MODULE thermcell + + Index: trunk/LMDZ.MARS/libf/phymars/thermcell_dqupdown.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/thermcell_dqupdown.F90 (revision 161) +++ trunk/LMDZ.MARS/libf/phymars/thermcell_dqupdown.F90 (revision 161) @@ -0,0 +1,301 @@ + subroutine thermcell_dqupdown(ngrid,nlay,ptimestep,fm0,entr0, & + & detr0,masse0,q_therm,dq_therm,ztvd,fm_down,ztv,charvar,lmax) + implicit none + +! #include "iniprint.h" +!======================================================================= +! +! Calcul du transport verticale dans la couche limite en presence +! de "thermiques" explicitement representes +! calcul du dq/dt une fois qu'on connait les ascendances +! Version modifiee pour prendre les downdrafts a la place de la +! subsidence compensatoire +!======================================================================= + +! ============================ INPUTS ============================ + + INTEGER, INTENT(IN) :: ngrid,nlay + REAL, INTENT(IN) :: ptimestep + REAL, INTENT(IN) :: fm0(ngrid,nlay+1) + REAL, INTENT(IN) :: entr0(ngrid,nlay),detr0(ngrid,nlay) + REAL, INTENT(IN) :: q_therm(ngrid,nlay) + REAL, INTENT(IN) :: fm_down(ngrid,nlay+1) + REAL, INTENT(IN) :: ztvd(ngrid,nlay) + REAL, INTENT(IN) :: ztv(ngrid,nlay) + CHARACTER (LEN=20), INTENT(IN) :: charvar + REAL, INTENT(IN) :: masse0(ngrid,nlay) + INTEGER, INTENT(IN) :: lmax(ngrid) + +! ============================ OUTPUTS =========================== + + REAL, INTENT(OUT) :: dq_therm(ngrid,nlay) + +! ============================ LOCAL ============================= + +! REAL detr0(ngrid,nlay) + REAL detrd(ngrid,nlay) + REAL entrd(ngrid,nlay) + REAL fmd(ngrid,nlay+1) + REAL q(ngrid,nlay) + REAL qa(ngrid,nlay) + REAL qd(ngrid,nlay) + INTEGER ig,k + LOGICAL active(ngrid,nlay) + INTEGER lmax_down(ngrid),lmin_down(ngrid) + INTEGER ncorec + +! =========== Init ============================================== + + entrd(:,:)=0. + detrd(:,:)=0. + qa(:,:)=q_therm(:,:) + q(:,:)=q_therm(:,:) + qd(:,:)=q_therm(:,:) + active(:,:)=.false. + +! previous calculation of zdthl_down uses the divergence of fmd +! so it can be negative without problem. Here we include the sign +! of fmd in the equations, so it has to be positive +! +! fmd(:,:)=-fm_down(:,:) +! +!! ========== Entrainment, Detrainement and Mass ================= +! +!! ========== DOWNDRAFT TRANSPORT DISABLED FOR NOW =============== +! +! do ig=1,ngrid +! if (ztv(ig,nlay)-ztvd(ig,nlay) .gt. 0.5) then +! print*,"downdraft non nul derniere couche !!! (dqupdown)" +! endif +! detrd(ig,nlay)=0. +! entrd(ig,nlay)=0. +! enddo +! +! do k=nlay-1,1,-1 +! do ig=1,ngrid +! +! if (ztv(ig,k)-ztvd(ig,k) .gt. 0.0001) then +! detrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztv(ig,k)-ztvd(ig,k+1))) & +! & /(ztv(ig,k)-ztvd(ig,k)) - fmd(ig,k)) +! entrd(ig,k)=MAX(0.,(fmd(ig,k+1)*(ztvd(ig,k)-ztvd(ig,k+1))) & +! & /(ztv(ig,k)-ztvd(ig,k))) +! +! endif +! enddo +! enddo +! +!! ======= We have computed entrainment and detrainment from a prescribed +!! mass flux and potential temp profile. Due to the way downdraft are parametrized, +!! this can yield negative entr and detr. We force it to be positive, but in +!! order to conserve tracers, we need to recompute an adequate mass flux +!! and modify interface rates, to preserve consistency. +! lmax_down(:)=1 +! lmin_down(:)=1 +! +! do k=1,nlay +! do ig=1,ngrid +! if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then +!! if (entrd(ig,k).gt.detrd(ig,k)) then +! lmax_down(ig)=min(k,lmax(ig)) +! endif +! enddo +! enddo +! do k=nlay,1,-1 +! do ig=1,ngrid +! if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then +!! if (detrd(ig,k).gt.entrd(ig,k)) then +! lmin_down(ig)=k +! endif +! enddo +! enddo +! +! +! fmd(:,:)=0. +! +! do ig=1,ngrid +! if ((lmax_down(ig).gt.1) .and. ((lmax_down(ig)-lmin_down(ig)).gt.1)) then +!! fmd(ig,lmax_down(ig))=0. +!! entrd(ig,lmax_down(ig))=detrd(ig,lmax_down(ig)) +!! detrd(ig,lmax_down(ig))=0. +!! print*,lmin_down(ig),lmax_down(ig),lmax(ig) +! +! fmd(ig,lmax_down(ig)+1)=0. +! +! do k=lmax_down(ig),lmin_down(ig)+1,-1 +! fmd(ig,k)=fmd(ig,k+1)+entrd(ig,k)-detrd(ig,k) +! enddo +! +! fmd(ig,lmin_down(ig))=0. +! detrd(ig,lmin_down(ig))=fmd(ig,lmin_down(ig)+1)+entrd(ig,lmin_down(ig)) +! +! else +! entrd(ig,:)=0. +! detrd(ig,:)=0. +! active(ig,:)=.false. +! endif +! +! enddo +! ncorec=0 +! do k=nlay,2,-1 +! do ig=1,ngrid +! if (fmd(ig,k).lt.0.) then +!! detrd(ig,k)=max(0.,detrd(ig,k)+fmd(ig,k-1)) +!! fmd(ig,k-1)=0. +!! entrd(ig,k-1)=0. +!! detrd(ig,k-1)=0. +!! lmin_down(ig)=k-1 +! fmd(ig,k)=fmd(ig,k+1) +! detrd(ig,k)=entrd(ig,k) +! ncorec=ncorec+1 +!! fmd(ig,k)=0. +!! detrd(ig,k)=entrd(ig,k)+fmd(ig,k+1) +! +! endif +! enddo +! enddo +! +! if (ncorec .ne. 0) then +! print*, 'corrections for negative downward mass flux :',ncorec +! endif +! print*, lmin_down(:),lmax_down(:) +! +! do k=2,nlay +! do ig=1,ngrid +! active(ig,k)=(k.ge.lmin_down(ig)).and.(k.le.lmax_down(ig)) & +! & .and.(((fmd(ig,k)+detrd(ig,k))*ptimestep).gt.1.e-6*masse0(ig,k)) +! enddo +! enddo +! +! +! do ig=1,ngrid +! do k=lmin_down(ig),lmax_down(ig) +! if(.not.active(ig,k)) then +! active(ig,:)=.false. +! endif +! enddo +! enddo +! +! if(charvar .eq. 'tke') then +! active(:,:)=.false. +! endif +! +!! do ig=1,ngrid +!! active(ig,lmax_down(ig))=(((fmd(ig,lmax_down(ig))+detrd(ig,lmax_down(ig)))* & +!! & ptimestep).gt.1.e-8*masse0(ig,lmax_down(ig))) +!! enddo +!! +!! do ig=1,ngrid +!! if (lmax_down(ig).gt.1) then +!! do k=lmax_down(ig)-1,lmin_down(ig),-1 +!! active(ig,k)=(((fmd(ig,k)+detrd(ig,k))* & +!! & ptimestep).gt.1.e-8*masse0(ig,k)) & +!! & .and. active(ig,k+1) +!! enddo +!! else +!! active(ig,:)=.false. +!! endif +!! enddo +!! +!! ========== qa : q in updraft ================================== +! + do k=2,nlay + do ig=1,ngrid + if ((fm0(ig,k+1)+detr0(ig,k))*ptimestep.gt. & + & 1.e-5*masse0(ig,k)) then + qa(ig,k)=(fm0(ig,k)*qa(ig,k-1)+entr0(ig,k)*q(ig,k)) & + & /(fm0(ig,k+1)+detr0(ig,k)) + + if ((qa(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then + print*,'qa<0!!!',charvar,ig,k,fm0(ig,k),qa(ig,k-1),entr0(ig,k),q(ig,k),fm0(ig,k+1),detr0(ig,k) + print*,'---------> Cancelling qa' + qa(ig,k)=q(ig,k) + endif + else + qa(ig,k)=q(ig,k) + endif + enddo + enddo + +! ========== qd : q in downdraft ================================= +! +! do k=nlay-1,1,-1 +! do ig=1,ngrid +! if (active(ig,k)) then +! qd(ig,k)=(fmd(ig,k+1)*qd(ig,k+1)+entrd(ig,k)*q(ig,k)) & +! & /(fmd(ig,k)+detrd(ig,k)) +! +! if ((qd(ig,k).lt.0.) .and. (charvar .ne. 'momentum')) then +! print*,'qd<0!!!',charvar,ig,k,fmd(ig,k),qd(ig,k),entrd(ig,k),q(ig,k),fmd(ig,k+1),detrd(ig,k),lmin_down(ig),lmax_down(ig) +! print*, '---------> cancelling qd, no downdraft for this gridpoint' +! qd(ig,k)=q(ig,k) +! active(ig,:)=.false. +! endif +! else +! qd(ig,k)=q(ig,k) +! endif +!! print*,'active,k,entr,detr,q,qd (down) :',active(ig,k),k,entrd(ig,k),detrd(ig,k),q(ig,k),qd(ig,k) +! enddo +! enddo +! +! +! ====== dq ====================================================== + + do ig=1,ngrid + if(active(ig,1)) then + + dq_therm(ig,1)=(detr0(ig,1)*qa(ig,1)+detrd(ig,1)*qd(ig,1) & + & +fm0(ig,2)*q(ig,2) & + & -entr0(ig,1)*q(ig,1)-entrd(ig,1)*q(ig,1) & + & -fmd(ig,2)*q(ig,1)) & + & *ptimestep/masse0(ig,1) + + else + dq_therm(ig,1)=(detr0(ig,1)*qa(ig,1)+fm0(ig,2)*q(ig,2) & + & -entr0(ig,1)*q(ig,1)) & + & *ptimestep/masse0(ig,1) + + endif + enddo + + do k=2,nlay-1 + do ig=1, ngrid + + if(active(ig,k)) then + + dq_therm(ig,k)=(detr0(ig,k)*qa(ig,k)+detrd(ig,k)*qd(ig,k) & + & +fm0(ig,k+1)*q(ig,k+1)+fmd(ig,k)*q(ig,k-1) & + & -entr0(ig,k)*q(ig,k)-entrd(ig,k)*q(ig,k) & + & -fm0(ig,k)*q(ig,k)-fmd(ig,k+1)*q(ig,k)) & + & *ptimestep/masse0(ig,k) + + + else + dq_therm(ig,k)=(detr0(ig,k)*qa(ig,k)+fm0(ig,k+1)*q(ig,k+1) & + & -entr0(ig,k)*q(ig,k)-fm0(ig,k)*q(ig,k)) & + & *ptimestep/masse0(ig,k) + + + endif + + enddo + enddo + + do ig=1, ngrid + + if(active(ig,nlay)) then + + dq_therm(ig,nlay)=(detr0(ig,nlay)*qa(ig,nlay)+detrd(ig,nlay)*qd(ig,nlay) & + & +fmd(ig,nlay)*q(ig,nlay-1) & + & -entr0(ig,nlay)*q(ig,nlay)-entrd(ig,nlay)*q(ig,nlay) & + & -fm0(ig,nlay)*q(ig,nlay)) & + & *ptimestep/masse0(ig,nlay) + + else + dq_therm(ig,nlay)=(detr0(ig,nlay)*qa(ig,nlay) & + & -entr0(ig,nlay)*q(ig,nlay)-fm0(ig,nlay)*q(ig,nlay)) & + & *ptimestep/masse0(ig,nlay) + endif + + enddo + return + end Index: trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90 (revision 161) +++ trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90 (revision 161) @@ -0,0 +1,1470 @@ +! +! + SUBROUTINE thermcell_main_mars(ngrid,nlay,nq,ptimestep & + & ,pplay,pplev,pphi,zlev,zlay & + & ,pu,pv,pt,pq,pq2 & + & ,pduadj,pdvadj,pdtadj,pdqadj,pdq2adj & + & ,fm,entr,detr,lmax & + & ,r_aspect & + & ,zw2,fraca & + & ,zpopsk,ztla,heatFlux,heatFlux_down & + & ,buoyancyOut, buoyancyEst) + + USE thermcell,ONLY:RG,RD,RKAPPA + IMPLICIT NONE + +!======================================================================= +! Mars version of thermcell_main. Author : A Colaitis +!======================================================================= +! ============== INPUTS ============== + + INTEGER, INTENT(IN) :: ngrid,nlay,nq + REAL, INTENT(IN) :: ptimestep,r_aspect + REAL, INTENT(IN) :: pt(ngrid,nlay) + REAL, INTENT(IN) :: pu(ngrid,nlay) + REAL, INTENT(IN) :: pv(ngrid,nlay) + REAL, INTENT(IN) :: pq(ngrid,nlay,nq) + REAL, INTENT(IN) :: pq2(ngrid,nlay) + REAL, INTENT(IN) :: pplay(ngrid,nlay) + REAL, INTENT(IN) :: pplev(ngrid,nlay+1) + REAL, INTENT(IN) :: pphi(ngrid,nlay) + REAL, INTENT(IN) :: zlay(ngrid,nlay) + REAL, INTENT(IN) :: zlev(ngrid,nlay+1) + +! ============== OUTPUTS ============== + + REAL, INTENT(OUT) :: pdtadj(ngrid,nlay) + REAL, INTENT(OUT) :: pduadj(ngrid,nlay) + REAL, INTENT(OUT) :: pdvadj(ngrid,nlay) + REAL, INTENT(OUT) :: pdqadj(ngrid,nlay,nq) + REAL, INTENT(OUT) :: pdq2adj(ngrid,nlay) + REAL, INTENT(OUT) :: zw2(ngrid,nlay+1) + +! Diagnostics + REAL, INTENT(OUT) :: heatFlux(ngrid,nlay) ! interface heatflux + REAL, INTENT(OUT) :: heatFlux_down(ngrid,nlay) ! interface heat flux from downdraft + REAL, INTENT(OUT) :: buoyancyOut(ngrid,nlay) ! interlayer buoyancy term + REAL, INTENT(OUT) :: buoyancyEst(ngrid,nlay) ! interlayer estimated buoyancy term + +! dummy variables when output not needed : + +! REAL :: heatFlux(ngrid,nlay) ! interface heatflux +! REAL :: heatFlux_down(ngrid,nlay) ! interface heat flux from downdraft +! REAL :: buoyancyOut(ngrid,nlay) ! interlayer buoyancy term +! REAL :: buoyancyEst(ngrid,nlay) ! interlayer estimated buoyancy term + + +! ============== LOCAL ================ + + INTEGER ig,k,l,ll,iq + INTEGER lmax(ngrid),lmin(ngrid),lalim(ngrid) + INTEGER lmix(ngrid) + INTEGER lmix_bis(ngrid) + REAL linter(ngrid) + REAL zmix(ngrid) + REAL zmax(ngrid) + REAL ztva(ngrid,nlay),zw_est(ngrid,nlay+1),ztva_est(ngrid,nlay) + REAL zmax_sec(ngrid) + REAL zh(ngrid,nlay) + REAL zdthladj(ngrid,nlay) + REAL zdthladj_down(ngrid,nlay) + REAL ztvd(ngrid,nlay) + REAL ztv(ngrid,nlay) + REAL zu(ngrid,nlay),zv(ngrid,nlay),zo(ngrid,nlay) + REAL zva(ngrid,nlay) + REAL zua(ngrid,nlay) + + REAL zta(ngrid,nlay) + REAL fraca(ngrid,nlay+1) + REAL q2(ngrid,nlay) + REAL rho(ngrid,nlay),rhobarz(ngrid,nlay),masse(ngrid,nlay) + REAL zpopsk(ngrid,nlay) + + REAL wmax(ngrid) + REAL wmax_sec(ngrid) + REAL fm(ngrid,nlay+1),entr(ngrid,nlay),detr(ngrid,nlay) + + REAL fm_down(ngrid,nlay+1) + + REAL ztla(ngrid,nlay) + + REAL f_star(ngrid,nlay+1),entr_star(ngrid,nlay) + REAL detr_star(ngrid,nlay) + REAL alim_star_tot(ngrid) + REAL alim_star(ngrid,nlay) + REAL alim_star_clos(ngrid,nlay) + REAL f(ngrid) + + REAL teta_th_int(ngrid,nlay) + REAL teta_env_int(ngrid,nlay) + REAL teta_down_int(ngrid,nlay) + + CHARACTER (LEN=20) :: modname + CHARACTER (LEN=80) :: abort_message + +! ============= PLUME VARIABLES ============ + + REAL w_est(ngrid,nlay+1) + REAL wa_moy(ngrid,nlay+1) + REAL wmaxa(ngrid) + REAL zdz,zbuoy(ngrid,nlay),zw2m + LOGICAL active(ngrid),activetmp(ngrid) + +! ========================================== + +! ============= HEIGHT VARIABLES =========== + + REAL num(ngrid) + REAL denom(ngrid) + REAL zlevinter(ngrid) + +! ========================================= + +! ============= DRY VARIABLES ============= + + REAL zw2_dry(ngrid,nlay+1) + REAL f_star_dry(ngrid,nlay+1) + REAL ztva_dry(ngrid,nlay+1) + REAL wmaxa_dry(ngrid) + REAL wa_moy_dry(ngrid,nlay+1) + REAL linter_dry(ngrid),zlevinter_dry(ngrid) + INTEGER lmix_dry(ngrid),lmax_dry(ngrid) + +! ========================================= + +! ============= CLOSURE VARIABLES ========= + + REAL zdenom(ngrid) + REAL alim_star2(ngrid) + REAL alim_star_tot_clos(ngrid) + INTEGER llmax + +! ========================================= + +! ============= FLUX2 VARIABLES =========== + + INTEGER ncorecfm1,ncorecfm2,ncorecfm3,ncorecalpha + INTEGER ncorecfm4,ncorecfm5,ncorecfm6,ncorecfm7,ncorecfm8 + REAL zfm + REAL f_old,ddd0,eee0,ddd,eee,zzz + REAL fomass_max,alphamax + +! ========================================= + +! ============= DTETA VARIABLES =========== + +! rien : on prends la divergence du flux turbulent + +! ========================================= + +! ============= DV2 VARIABLES ============= +! not used for now + + REAL qa(ngrid,nlay),detr_dv2(ngrid,nlay),zf,zf2 + REAL wvd(ngrid,nlay+1),wud(ngrid,nlay+1) + REAL gamma0(ngrid,nlay+1),gamma(ngrid,nlay+1) + REAL ue(ngrid,nlay),ve(ngrid,nlay) + LOGICAL ltherm(ngrid,nlay) + REAL dua(ngrid,nlay),dva(ngrid,nlay) + INTEGER iter + INTEGER nlarga0 + +! ========================================= + +!----------------------------------------------------------------------- +! initialisation: +! --------------- + + zu(:,:)=pu(:,:) + zv(:,:)=pv(:,:) + ztv(:,:)=pt(:,:)/zpopsk(:,:) + +!------------------------------------------------------------------------ +! -------------------- +! +! +! + + + + + + + + + + + +! +! +! wa, fraca, wd, fracd -------------------- zlev(2), rhobarz +! wh,wt,wo ... +! +! + + + + + + + + + + + zh,zu,zv,zo,rho +! +! +! -------------------- zlev(1) +! \\\\\\\\\\\\\\\\\\\\ +! +! + +!----------------------------------------------------------------------- +! Calcul des altitudes des couches +!----------------------------------------------------------------------- + +! do l=2,nlay +! zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG +! enddo +! zlev(:,1)=0. +! zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG + +! zlay(:,:)=pphi(:,:)/RG +!----------------------------------------------------------------------- +! Calcul des densites +!----------------------------------------------------------------------- + + rho(:,:)=pplay(:,:)/(RD*pt(:,:)) + + rhobarz(:,1)=rho(:,1) + + do l=2,nlay + rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1)) + enddo + +!calcul de la masse + do l=1,nlay + masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG + enddo + + +!------------------------------------------------------------------ +! +! /|\ +! -------- | F_k+1 ------- +! ----> D_k +! /|\ <---- E_k , A_k +! -------- | F_k --------- +! ----> D_k-1 +! <---- E_k-1 , A_k-1 +! +! +! --------------------------- +! +! ----- F_lmax+1=0 ---------- \ +! lmax (zmax) | +! --------------------------- | +! | +! --------------------------- | +! | +! --------------------------- | +! | +! --------------------------- | +! | +! --------------------------- | +! | E +! --------------------------- | D +! | +! --------------------------- | +! | +! --------------------------- \ | +! lalim | | +! --------------------------- | | +! | | +! --------------------------- | | +! | A | +! --------------------------- | | +! | | +! --------------------------- | | +! lmin (=1 pour le moment) | | +! ----- F_lmin=0 ------------ / / +! +! --------------------------- +! ////////////////////////// +! + +!============================================================================= +! Calculs initiaux ne faisant pas intervenir les changements de phase +!============================================================================= + +!------------------------------------------------------------------ +! 1. alim_star est le profil vertical de l'alimentation a la base du +! panache thermique, calcule a partir de la flotabilite de l'air sec +! 2. lmin et lalim sont les indices inferieurs et superieurs de alim_star +!------------------------------------------------------------------ +! + entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0. + lmin=1 + +!----------------------------------------------------------------------------- +! 3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un +! panache sec conservatif (e=d=0) alimente selon alim_star +! Il s'agit d'un calcul de type CAPE +! zmax_sec est utilise pour determiner la geometrie du thermique. +!------------------------------------------------------------------------------ +!--------------------------------------------------------------------------------- +!calcul du melange et des variables dans le thermique +!-------------------------------------------------------------------------------- + +! =========================================================================== +! ===================== PLUME =============================================== +! =========================================================================== + +! Initialisations des variables reeles + ztva(:,:)=ztv(:,:) + ztva_est(:,:)=ztva(:,:) + ztla(:,:)=0. + + zbuoy(:,:)=0. + w_est(:,:)=0. + f_star(:,:)=0. + wa_moy(:,:)=0. + linter(:)=1. +! Initialisation des variables entieres + lmix(:)=1 + lmix_bis(:)=2 + wmaxa(:)=0. + lalim(:)=1 + +!------------------------------------------------------------------------- +! On ne considere comme actif que les colonnes dont les deux premieres +! couches sont instables. +!------------------------------------------------------------------------- + active(:)=ztv(:,1)>ztv(:,2) + do ig=1,ngrid + if (ztv(ig,1)>=(ztv(ig,2))) then + alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.) & + & *sqrt(zlev(ig,2)) + lalim(ig)=2 + alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,1) + endif + enddo + + do l=2,nlay-1 + do ig=1,ngrid + if (ztv(ig,l)>(ztv(ig,l+1)+0.5) .and. ztv(ig,1)>=ztv(ig,l) .and. ztv(ig,l-1)>(ztv(ig,l))) then + alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.) & + & *sqrt(zlev(ig,l+1)) + lalim(ig)=l+1 + alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) + endif + enddo + enddo + do l=1,nlay + do ig=1,ngrid + if (alim_star_tot(ig) > 1.e-10 ) then + alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig) + endif + enddo + enddo + + alim_star_tot(:)=1. + +!------------------------------------------------------------------------------ +! Calcul dans la premiere couche +! On decide dans cette version que le thermique n'est actif que si la premiere +! couche est instable. +! Pourrait etre change si on veut que le thermiques puisse se déclencher +! dans une couche l>1 +!------------------------------------------------------------------------------ + + do ig=1,ngrid +! Le panache va prendre au debut les caracteristiques de l'air contenu +! dans cette couche. + if (active(ig)) then + ztla(ig,1)=ztv(ig,1) +!cr: attention, prise en compte de f*(1)=1 => AC : what ? f*(1) =0. ! (d'ou f*(2)=a*(1) +! dans un panache conservatif + f_star(ig,1)=0. + f_star(ig,2)=alim_star(ig,1) + zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2) & + & *(zlev(ig,2)-zlev(ig,1)) & + & *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) + w_est(ig,2)=zw2(ig,2) + + endif + enddo + + +!============================================================================== +!boucle de calcul de la vitesse verticale dans le thermique +!============================================================================== + do l=2,nlay-1 +!============================================================================== + + +! On decide si le thermique est encore actif ou non +! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test + do ig=1,ngrid + active(ig)=active(ig) & + & .and. zw2(ig,l)>1.e-10 & + & .and. f_star(ig,l)+alim_star(ig,l)>1.e-10 + enddo + +!--------------------------------------------------------------------------- +! calcul des proprietes thermodynamiques et de la vitesse de la couche l +! sans tenir compte du detrainement et de l'entrainement dans cette +! couche +! C'est a dire qu'on suppose +! ztla(l)=ztla(l-1) +! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer +! avant) a l'alimentation pour avoir un calcul plus propre +!--------------------------------------------------------------------------- + + do ig=1,ngrid + if(active(ig)) then + + if(l .lt. lalim(ig)) then + ztva_est(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & + & alim_star(ig,l)*ztv(ig,l)) & + & /(f_star(ig,l)+alim_star(ig,l)) + else + ztva_est(ig,l)=ztla(ig,l-1) + endif + + zdz=zlev(ig,l+1)-zlev(ig,l) + zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) + if (((2.5*zbuoy(ig,l)/w_est(ig,l)-0.0015) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then + w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*0.0015 & + & -2.*zdz*w_est(ig,l)*0.045*(2.5*zbuoy(ig,l)/w_est(ig,l)-0.0015)**0.6) + else + w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*0.0015) + endif + if (w_est(ig,l+1).lt.0.) then + w_est(ig,l+1)=zw2(ig,l) + endif + endif + enddo + +!------------------------------------------------- +!calcul des taux d'entrainement et de detrainement +!------------------------------------------------- + + do ig=1,ngrid + if (active(ig)) then + + zw2m=w_est(ig,l+1) + if((2.5*(zbuoy(ig,l)/zw2m)-0.0015) .gt. 0.) then + entr_star(ig,l)=f_star(ig,l)*zdz* & +! & 0.0118*((zbuoy(ig,l)/zw2m)/0.043)**(1./1.65) +! & 0.0090*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.90) +! & 0.0120*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.6) + & MAX(0.,0.045*(2.5*(zbuoy(ig,l)/zw2m)-0.0015)**0.6) + else + entr_star(ig,l)=0. + endif + if(zbuoy(ig,l) .gt. 0.) then + if(l .lt. lalim(ig)) then + detr_star(ig,l)=0. + else +! detr_star(ig,l) = f_star(ig,l)*zdz* & +! & 0.0105*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.7) +! detr_star(ig,l) = f_star(ig,l)*zdz* & +! & 0.0085*((zbuoy(ig,l)/zw2m)/0.05)**(1./1.55) + +! last baseline from direct les +! detr_star(ig,l) = f_star(ig,l)*zdz* & +! & 0.065*(2.5*(zbuoy(ig,l)/zw2m))**0.75 + +! new param from continuity eq with a fit on dfdz + detr_star(ig,l) = f_star(ig,l)*zdz* & + & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) + +! & 0.014*((zbuoy(ig,l)/zw2m)/0.05)**(1./1.35) +! detr_star(ig,l) = f_star(ig,l)*zdz* & +! & ((zbuoy(ig,l)/zw2m)/2.222)! + 0.0002) + + endif + else + detr_star(ig,l)=f_star(ig,l)*zdz* & + & MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005) + +! & *5.*(-afact*zbetalpha*zbuoy(ig,l)/zw2m) +! & *5.*(-afact*zbuoy(ig,l)/zw2m) + +! last baseline from direct les +! & 0.065*(-2.5*(zbuoy(ig,l)/zw2m))**0.75 + +! new param from continuity eq with a fit on dfdz + + + endif + +! En dessous de lalim, on prend le max de alim_star et entr_star pour +! alim_star et 0 sinon + + if (l.lt.lalim(ig)) then + alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l)) + entr_star(ig,l)=0. + endif + +! Calcul du flux montant normalise + + f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l) & + & -detr_star(ig,l) + + endif + enddo + + +!---------------------------------------------------------------------------- +!calcul de la vitesse verticale en melangeant Tl et qt du thermique +!--------------------------------------------------------------------------- + + activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10 + do ig=1,ngrid + if (activetmp(ig)) then + + ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+ & + & (alim_star(ig,l)+entr_star(ig,l))*ztv(ig,l)) & + & /(f_star(ig,l+1)+detr_star(ig,l)) + + endif + enddo + + do ig=1,ngrid + if (activetmp(ig)) then + ztva(ig,l) = ztla(ig,l) + zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) + + if (((2.5*zbuoy(ig,l)/zw2(ig,l)-0.0015) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then + zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*2.5*zbuoy(ig,l)- & + & 2.*zdz*zw2(ig,l)*0.0015-2.*zdz*zw2(ig,l)*0.045*(2.5*zbuoy(ig,l)/zw2(ig,l)-0.0015)**0.6) + else + zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*0.0015) + endif + endif + enddo + +!--------------------------------------------------------------------------- +!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max +!--------------------------------------------------------------------------- + + do ig=1,ngrid + if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then + print*,'On tombe sur le cas particulier de thermcell_plume' + zw2(ig,l+1)=0. + linter(ig)=l+1 + endif + + if (zw2(ig,l+1).lt.0.) then + linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l)) & + & -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l)) + zw2(ig,l+1)=0. + endif + wa_moy(ig,l+1)=sqrt(zw2(ig,l+1)) + + if (wa_moy(ig,l+1).gt.wmaxa(ig)) then +! lmix est le niveau de la couche ou w (wa_moy) est maximum +!on rajoute le calcul de lmix_bis + lmix(ig)=l+1 + wmaxa(ig)=wa_moy(ig,l+1) + endif + enddo + +!========================================================================= +! FIN DE LA BOUCLE VERTICALE + enddo +!========================================================================= + +!on recalcule alim_star_tot + do ig=1,ngrid + alim_star_tot(ig)=0. + enddo + do ig=1,ngrid + do l=1,lalim(ig)-1 + alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l) + enddo + enddo + + +! =========================================================================== +! ================= FIN PLUME =============================================== +! =========================================================================== + + +! =========================================================================== +! ================= HEIGHT ================================================== +! =========================================================================== + +! Attention, w2 est transforme en sa racine carree dans cette routine + +!------------------------------------------------------------------------------- +! Calcul des caracteristiques du thermique:zmax,zmix,wmax +!------------------------------------------------------------------------------- + +!calcul de la hauteur max du thermique + do ig=1,ngrid + lmax(ig)=lalim(ig) + enddo + do ig=1,ngrid + do l=nlay,lalim(ig)+1,-1 + if (zw2(ig,l).le.1.e-10) then + lmax(ig)=l-1 + endif + enddo + enddo + +! On traite le cas particulier qu'il faudrait éviter ou le thermique +! atteind le haut du modele ... + do ig=1,ngrid + if ( zw2(ig,nlay) > 1.e-10 ) then + print*,'WARNING !!!!! W2 thermiques non nul derniere couche ' + lmax(ig)=nlay + endif + enddo + +! pas de thermique si couche 1 stable +! do ig=1,ngrid +! if (lmin(ig).gt.1) then +! lmax(ig)=1 +! lmin(ig)=1 +! lalim(ig)=1 +! endif +! enddo +! +! Determination de zw2 max + do ig=1,ngrid + wmax(ig)=0. + enddo + + do l=1,nlay + do ig=1,ngrid + if (l.le.lmax(ig)) then + if (zw2(ig,l).lt.0.)then + print*,'pb2 zw2<0' + endif + zw2(ig,l)=sqrt(zw2(ig,l)) + wmax(ig)=max(wmax(ig),zw2(ig,l)) + else + zw2(ig,l)=0. + endif + enddo + enddo +! Longueur caracteristique correspondant a la hauteur des thermiques. + do ig=1,ngrid + zmax(ig)=0. + zlevinter(ig)=zlev(ig,1) + enddo + + num(:)=0. + denom(:)=0. + do ig=1,ngrid + do l=1,nlay + num(ig)=num(ig)+zw2(ig,l)*zlev(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) + denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) + enddo + enddo + do ig=1,ngrid + if (denom(ig).gt.1.e-10) then + zmax(ig)=2.*num(ig)/denom(ig) + endif + enddo + +! def de zmix continu (profil parabolique des vitesses) + do ig=1,ngrid + if (lmix(ig).gt.1) then +! test + if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) & + & *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) & + & -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) & + & *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10) & + & then +! + zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) & + & *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2) & + & -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) & + & *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2)) & + & /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig))) & + & *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1))) & + & -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1)) & + & *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig)))))) + else + zmix(ig)=zlev(ig,lmix(ig)) + print*,'pb zmix' + endif + else + zmix(ig)=0. + endif +!test + if ((zmax(ig)-zmix(ig)).le.0.) then + zmix(ig)=0.9*zmax(ig) + endif + enddo +! +! calcul du nouveau lmix correspondant + do ig=1,ngrid + do l=1,nlay + if (zmix(ig).ge.zlev(ig,l).and. & + & zmix(ig).lt.zlev(ig,l+1)) then + lmix(ig)=l + endif + enddo + enddo + + +! Attention, w2 est transforme en sa racine carree dans cette routine + +! =========================================================================== +! ================= FIN HEIGHT ============================================== +! =========================================================================== + + +! =========================================================================== +! ================= DRY ===================================================== +! =========================================================================== + +!initialisations + do ig=1,ngrid + do l=1,nlay+1 + zw2_dry(ig,l)=0. + wa_moy_dry(ig,l)=0. + enddo + enddo + do ig=1,ngrid + do l=1,nlay + ztva_dry(ig,l)=ztv(ig,l) + enddo + enddo + do ig=1,ngrid + wmax_sec(ig)=0. + wmaxa_dry(ig)=0. + enddo +!calcul de la vitesse a partir de la CAPE en melangeant thetav + + +! Calcul des F^*, integrale verticale de E^* + f_star_dry(:,1)=0. + do l=1,nlay + f_star_dry(:,l+1)=f_star_dry(:,l)+alim_star(:,l) + enddo + +! niveau (reel) auquel zw2 s'annule FH :n'etait pas initialise + linter_dry(:)=0. + +! couche la plus haute concernee par le thermique. + lmax_dry(:)=1 + +! Le niveau linter est une variable continue qui se trouve dans la couche +! lmax + + do l=1,nlay-2 + do ig=1,ngrid + if (l.eq.lmin(ig).and.lalim(ig).gt.1) then + +!------------------------------------------------------------------------ +! Calcul de la vitesse en haut de la premiere couche instable. +! Premiere couche du panache thermique +!------------------------------------------------------------------------ + + zw2_dry(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1) & + & *(zlev(ig,l+1)-zlev(ig,l)) & + & *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l)) +!------------------------------------------------------------------------ +! Tant que la vitesse en bas de la couche et la somme du flux de masse +! et de l'entrainement (c'est a dire le flux de masse en haut) sont +! positifs, on calcul +! 1. le flux de masse en haut f_star(ig,l+1) +! 2. la temperature potentielle virtuelle dans la couche ztva(ig,l) +! 3. la vitesse au carré en haut zw2(ig,l+1) +!------------------------------------------------------------------------ + + else if (zw2_dry(ig,l).ge.1e-10) then + + ztva_dry(ig,l)=(f_star_dry(ig,l)*ztva_dry(ig,l-1)+alim_star(ig,l) & + & *ztv(ig,l))/f_star_dry(ig,l+1) + zw2_dry(ig,l+1)=zw2_dry(ig,l)*(f_star_dry(ig,l)/f_star_dry(ig,l+1))**2+ & + & 2.*RG*(ztva_dry(ig,l)-ztv(ig,l))/ztv(ig,l) & + & *(zlev(ig,l+1)-zlev(ig,l)) + endif +! determination de zmax continu par interpolation lineaire +!------------------------------------------------------------------------ + + if (zw2_dry(ig,l+1)>0. .and. zw2_dry(ig,l+1).lt.1.e-10) then +! stop'On tombe sur le cas particulier de thermcell_dry' +! print*,'On tombe sur le cas particulier de thermcell_dry' + zw2_dry(ig,l+1)=0. + linter_dry(ig)=l+1 + lmax_dry(ig)=l + endif + + if (zw2_dry(ig,l+1).lt.0.) then + linter_dry(ig)=(l*(zw2_dry(ig,l+1)-zw2_dry(ig,l)) & + & -zw2_dry(ig,l))/(zw2_dry(ig,l+1)-zw2_dry(ig,l)) + zw2_dry(ig,l+1)=0. + lmax_dry(ig)=l + endif + + wa_moy_dry(ig,l+1)=sqrt(zw2_dry(ig,l+1)) + + if (wa_moy_dry(ig,l+1).gt.wmaxa_dry(ig)) then +! lmix est le niveau de la couche ou w (wa_moy) est maximum + lmix_dry(ig)=l+1 + wmaxa_dry(ig)=wa_moy_dry(ig,l+1) + endif + enddo + enddo +! +! Determination de zw2 max + do ig=1,ngrid + wmax_sec(ig)=0. + enddo + do l=1,nlay + do ig=1,ngrid + if (l.le.lmax_dry(ig)) then + zw2_dry(ig,l)=sqrt(zw2_dry(ig,l)) + wmax_sec(ig)=max(wmax_sec(ig),zw2_dry(ig,l)) + else + zw2_dry(ig,l)=0. + endif + enddo + enddo + +! Longueur caracteristique correspondant a la hauteur des thermiques. + do ig=1,ngrid + zmax_sec(ig)=0. + zlevinter_dry(ig)=zlev(ig,1) + enddo + do ig=1,ngrid +! calcul de zlevinter + zlevinter_dry(ig)=zlev(ig,lmax_dry(ig)) + & + & (linter_dry(ig)-lmax_dry(ig))*(zlev(ig,lmax_dry(ig)+1)-zlev(ig,lmax_dry(ig))) + zmax_sec(ig)=max(zmax_sec(ig),zlevinter_dry(ig)-zlev(ig,lmin(ig))) + enddo + +! =========================================================================== +! ============= FIN DRY ===================================================== +! =========================================================================== + + +! Choix de la fonction d'alimentation utilisee pour la fermeture. + + alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:) + +! =========================================================================== +! ============= CLOSURE ===================================================== +! =========================================================================== + +!------------------------------------------------------------------------------- +! Fermeture,determination de f +!------------------------------------------------------------------------------- +! Appel avec la version seche + + alim_star2(:)=0. + alim_star_tot_clos(:)=0. + f(:)=0. + +! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine + llmax=1 + do ig=1,ngrid + if (lalim(ig)>llmax) llmax=lalim(ig) + enddo + + +! Calcul des integrales sur la verticale de alim_star et de +! alim_star^2/(rho dz) + do k=1,llmax-1 + do ig=1,ngrid + if (k1.e-10) then + f(ig)=wmax_sec(ig)*alim_star_tot_clos(ig)/ & + & (max(500.,zmax_sec(ig))*r_aspect*alim_star2(ig)) + endif + enddo + +! =========================================================================== +! ============= FIN CLOSURE ================================================= +! =========================================================================== + + +! =========================================================================== +! ============= FLUX2 ======================================================= +! =========================================================================== + +!------------------------------------------------------------------------------- +!deduction des flux +!------------------------------------------------------------------------------- + + fomass_max=0.8 + alphamax=0.5 + + ncorecfm1=0 + ncorecfm2=0 + ncorecfm3=0 + ncorecfm4=0 + ncorecfm5=0 + ncorecfm6=0 + ncorecfm7=0 + ncorecfm8=0 + ncorecalpha=0 + +!------------------------------------------------------------------------- +! Multiplication par le flux de masse issu de la femreture +!------------------------------------------------------------------------- + + do l=1,nlay + entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l)) + detr(:,l)=f(:)*detr_star(:,l) + enddo + + do l=1,nlay + do ig=1,ngrid + if (l.lt.lmax(ig)) then + fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l) + elseif(l.eq.lmax(ig)) then + fm(ig,l+1)=0. + detr(ig,l)=fm(ig,l)+entr(ig,l) + else + fm(ig,l+1)=0. + endif + enddo + enddo + +! Test provisoire : pour comprendre pourquoi on corrige plein de fois +! le cas fm6, on commence par regarder une premiere fois avant les +! autres corrections. + + do l=1,nlay + do ig=1,ngrid + if (detr(ig,l).gt.fm(ig,l)) then + ncorecfm8=ncorecfm8+1 + endif + enddo + enddo + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! FH Version en cours de test; +! par rapport a thermcell_flux, on fait une grande boucle sur "l" +! et on modifie le flux avec tous les contr�les appliques d'affilee +! pour la meme couche +! Momentanement, on duplique le calcule du flux pour pouvoir comparer +! les flux avant et apres modif +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + do l=1,nlay + + do ig=1,ngrid + if (l.lt.lmax(ig)) then + fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l) + elseif(l.eq.lmax(ig)) then + fm(ig,l+1)=0. + detr(ig,l)=fm(ig,l)+entr(ig,l) + else + fm(ig,l+1)=0. + endif + enddo + + +!------------------------------------------------------------------------- +! Verification de la positivite des flux de masse +!------------------------------------------------------------------------- + + do ig=1,ngrid + if (fm(ig,l+1).lt.0.) then + print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1) + ncorecfm1=ncorecfm1+1 + fm(ig,l+1)=fm(ig,l) + detr(ig,l)=entr(ig,l) + endif + enddo + +! Les "optimisations" du flux sont desactivees : moins de bidouilles +! je considere que celles ci ne sont pas justifiees ou trop delicates +! pour MARS, d'apres des observations LES. +!------------------------------------------------------------------------- +!Test sur fraca croissant +!------------------------------------------------------------------------- +! if (iflag_thermals_optflux==0) then +! do ig=1,ngrid +! if (l.ge.lalim(ig).and.l.le.lmax(ig) & +! & .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then +!! zzz est le flux en l+1 a frac constant +! zzz=fm(ig,l)*rhobarz(ig,l+1)*zw2(ig,l+1) & +! & /(rhobarz(ig,l)*zw2(ig,l)) +! if (fm(ig,l+1).gt.zzz) then +! detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz +! fm(ig,l+1)=zzz +! ncorecfm4=ncorecfm4+1 +! endif +! endif +! enddo +! endif +! +!------------------------------------------------------------------------- +!test sur flux de masse croissant +!------------------------------------------------------------------------- +! if (iflag_thermals_optflux==0) then +! do ig=1,ngrid +! if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then +! f_old=fm(ig,l+1) +! fm(ig,l+1)=fm(ig,l) +! detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1) +! ncorecfm5=ncorecfm5+1 +! endif +! enddo +! endif +! +!------------------------------------------------------------------------- +!detr ne peut pas etre superieur a fm +!------------------------------------------------------------------------- + + do ig=1,ngrid + if (detr(ig,l).gt.fm(ig,l)) then + ncorecfm6=ncorecfm6+1 + detr(ig,l)=fm(ig,l) + entr(ig,l)=fm(ig,l+1) + +! Dans le cas ou on est au dessus de la couche d'alimentation et que le +! detrainement est plus fort que le flux de masse, on stope le thermique. + endif + + if(l.gt.lmax(ig)) then + detr(ig,l)=0. + fm(ig,l+1)=0. + entr(ig,l)=0. + endif + enddo + +!------------------------------------------------------------------------- +!fm ne peut pas etre negatif +!------------------------------------------------------------------------- + + do ig=1,ngrid + if (fm(ig,l+1).lt.0.) then + detr(ig,l)=detr(ig,l)+fm(ig,l+1) + fm(ig,l+1)=0. + ncorecfm2=ncorecfm2+1 + endif + enddo + +!----------------------------------------------------------------------- +!la fraction couverte ne peut pas etre superieure a 1 +!----------------------------------------------------------------------- +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! FH Partie a revisiter. +! Il semble qu'etaient codees ici deux optiques dans le cas +! F/ (rho *w) > 1 +! soit limiter la hauteur du thermique en considerant que c'est +! la derniere chouche, soit limiter F a rho w. +! Dans le second cas, il faut en fait limiter a un peu moins +! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin +! dans thermcell_main et qu'il semble de toutes facons deraisonable +! d'avoir des fractions de 1.. +! Ci dessous, et dans l'etat actuel, le premier des deux if est +! sans doute inutile. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + do ig=1,ngrid + if (zw2(ig,l+1).gt.1.e-10) then + zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax + if ( fm(ig,l+1) .gt. zfm) then + f_old=fm(ig,l+1) + fm(ig,l+1)=zfm + detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1) + ncorecalpha=ncorecalpha+1 + endif + endif + + enddo + +! Fin de la grande boucle sur les niveaux verticaux + enddo + +!----------------------------------------------------------------------- +! On fait en sorte que la quantite totale d'air entraine dans le +! panache ne soit pas trop grande comparee a la masse de la maille +!----------------------------------------------------------------------- + + do l=1,nlay-1 + do ig=1,ngrid + eee0=entr(ig,l) + ddd0=detr(ig,l) + eee=entr(ig,l)-masse(ig,l)*fomass_max/ptimestep + ddd=detr(ig,l)-eee + if (eee.gt.0.) then + ncorecfm3=ncorecfm3+1 + entr(ig,l)=entr(ig,l)-eee + if ( ddd.gt.0.) then +! l'entrainement est trop fort mais l'exces peut etre compense par une +! diminution du detrainement) + detr(ig,l)=ddd + else +! l'entrainement est trop fort mais l'exces doit etre compense en partie +! par un entrainement plus fort dans la couche superieure + if(l.eq.lmax(ig)) then + detr(ig,l)=fm(ig,l)+entr(ig,l) + else + entr(ig,l+1)=entr(ig,l+1)-ddd + detr(ig,l)=0. + fm(ig,l+1)=fm(ig,l)+entr(ig,l) + detr(ig,l)=0. + endif + endif + endif + enddo + enddo +! +! ddd=detr(ig)-entre +!on s assure que tout s annule bien en zmax + do ig=1,ngrid + fm(ig,lmax(ig)+1)=0. + entr(ig,lmax(ig))=0. + detr(ig,lmax(ig))=fm(ig,lmax(ig))+entr(ig,lmax(ig)) + enddo + +!----------------------------------------------------------------------- +! Impression du nombre de bidouilles qui ont ete necessaires +!----------------------------------------------------------------------- + +!IM 090508 beg + if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > ngrid/4. ) then + print*,'thermcell warning : large number of corrections' + print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',& + & ncorecfm2,'x fm2',ncorecfm3,'x fm3 et', & + & ncorecfm4,'x fm4',ncorecfm5,'x fm5 et', & + & ncorecfm6,'x fm6', & + & ncorecfm7,'x fm7', & + & ncorecfm8,'x fm8', & + & ncorecalpha,'x alpha' + endif + +! =========================================================================== +! ============= FIN FLUX2 =================================================== +! =========================================================================== + + +! =========================================================================== +! ============= TRANSPORT =================================================== +! =========================================================================== + +!------------------------------------------------------------------ +! calcul du transport vertical +!------------------------------------------------------------------ + +! ------------------------------------------------------------------ +! Transport de teta dans l'updraft : (remplace thermcell_dq) +! ------------------------------------------------------------------ + + zdthladj(:,:)=0. + + do ig=1,ngrid + if(lmax(ig) .gt. 1) then + do k=1,lmax(ig) + zdthladj(ig,k)=(1./masse(ig,k))*(fm(ig,k+1)*ztv(ig,k+1)- & + & fm(ig,k)*ztv(ig,k)+fm(ig,k)*ztva(ig,k)-fm(ig,k+1)*ztva(ig,k+1)) + if (ztv(ig,k) + ptimestep*zdthladj(ig,k) .le. 0.) then + print*,'q<0 in thermcell_dTeta up: qenv .. dq : ', ztv(ig,k),ptimestep*zdthladj(ig,k) + endif + enddo + endif + enddo + +! ------------------------------------------------------------------ +! Prescription des proprietes du downdraft +! ------------------------------------------------------------------ + + ztvd(:,:)=ztv(:,:) + fm_down(:,:)=0. + do ig=1,ngrid + if (lmax(ig) .gt. 1) then + do l=1,lmax(ig) + if(zlay(ig,l) .le. 0.7*zmax(ig)) then + fm_down(ig,l) =-1.8*fm(ig,l) + endif + + if(zlay(ig,l) .le. 0.06*zmax(ig)) then + ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(sqrt((zlay(ig,l)/zmax(ig))/0.122449) - 1.)*(ztva(ig,l)/ztv(ig,l) - 1.))) + elseif(zlay(ig,l) .le. 0.6*zmax(ig)) then + ztvd(ig,l)=ztv(ig,l)*max(0.,1.-0.3*(ztva(ig,l)/ztv(ig,l) - 1.)) + elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then + ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/0.333333)*(ztva(ig,l)/ztv(ig,l) - 1.))) + else + ztvd(ig,l)=ztv(ig,l) + endif + +! if(zlay(ig,l) .le. 0.06*zmax(ig)) then +! ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig))/19.231 + 0.9938) +! elseif(zlay(ig,l) .le. 0.6*zmax(ig)) then +! ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig)-0.075)/187.931 + 0.9982) +! else +!! elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then +! ztvd(ig,l)=ztv(ig,l)*max(0.,(zlay(ig,l)/zmax(ig) -0.60)/(-1333) + 1.00025) +!! else +!! ztvd(ig,l)=ztv(ig,l) +! endif + + enddo + endif + enddo + +! ------------------------------------------------------------------ +! Transport de teta dans le downdraft : (remplace thermcell_dq) +! ------------------------------------------------------------------ + + zdthladj_down(:,:)=0. + + do ig=1,ngrid + if(lmax(ig) .gt. 1) then + do k=1,lmax(ig) + zdthladj_down(ig,k)=(1./masse(ig,k))*(fm_down(ig,k+1)*ztv(ig,k+1)- & + & fm_down(ig,k)*ztv(ig,k)+fm_down(ig,k)*ztvd(ig,k)-fm_down(ig,k+1)*ztvd(ig,k+1)) + if (ztv(ig,k) + ptimestep*zdthladj_down(ig,k) .le. 0.) then + print*,'q<0 in thermcell_dTeta down: qenv .. dq : ', ztv(ig,k),ptimestep*zdthladj_down(ig,k) + endif + enddo + endif + enddo + +!------------------------------------------------------------------ +! Calcul de la fraction de l'ascendance +!------------------------------------------------------------------ + do ig=1,ngrid + fraca(ig,1)=0. + fraca(ig,nlay+1)=0. + enddo + do l=2,nlay + do ig=1,ngrid + if (zw2(ig,l).gt.1.e-10) then + fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l)) + else + fraca(ig,l)=0. + endif + enddo + enddo + + + +! =========================================================================== +! ============= DV2 ========================================================= +! =========================================================================== +! ROUTINE OVERIDE : ne prends pas en compte le downdraft +! de plus, le gradient de pression horizontal semble tout deregler... A VOIR + + if (0 .eq. 1) then + +!------------------------------------------------------------------ +! calcul du transport vertical du moment horizontal +!------------------------------------------------------------------ + +! Calcul du transport de V tenant compte d'echange par gradient +! de pression horizontal avec l'environnement + +! calcul du detrainement +!--------------------------- + + nlarga0=0. + + do k=1,nlay + do ig=1,ngrid + detr_dv2(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k) + enddo + enddo + +! calcul de la valeur dans les ascendances + do ig=1,ngrid + zua(ig,1)=zu(ig,1) + zva(ig,1)=zv(ig,1) + ue(ig,1)=zu(ig,1) + ve(ig,1)=zv(ig,1) + enddo + + gamma(1:ngrid,1)=0. + do k=2,nlay + do ig=1,ngrid + ltherm(ig,k)=(fm(ig,k+1)+detr_dv2(ig,k))*ptimestep > 1.e-5*masse(ig,k) + if(ltherm(ig,k).and.zmax(ig)>0.) then + gamma0(ig,k)=masse(ig,k) & + & *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) ) & + & *0.5/zmax(ig) & + & *1. + else + gamma0(ig,k)=0. + endif + if (ltherm(ig,k).and.zmax(ig)<=0.) nlarga0=nlarga0+1 + enddo + enddo + + gamma(:,:)=0. + + do k=2,nlay + + do ig=1,ngrid + + if (ltherm(ig,k)) then + dua(ig,k)=zua(ig,k-1)-zu(ig,k-1) + dva(ig,k)=zva(ig,k-1)-zv(ig,k-1) + else + zua(ig,k)=zu(ig,k) + zva(ig,k)=zv(ig,k) + ue(ig,k)=zu(ig,k) + ve(ig,k)=zv(ig,k) + endif + enddo + + +! Debut des iterations +!---------------------- + +! AC WARNING : SALE ! + + do iter=1,5 + do ig=1,ngrid +! Pour memoire : calcul prenant en compte la fraction reelle +! zf=0.5*(fraca(ig,k)+fraca(ig,k+1)) +! zf2=1./(1.-zf) +! Calcul avec fraction infiniement petite + zf=0. + zf2=1. + +! la première fois on multiplie le coefficient de freinage +! par le module du vent dans la couche en dessous. +! Mais pourquoi donc ??? + if (ltherm(ig,k)) then +! On choisit une relaxation lineaire. +! gamma(ig,k)=gamma0(ig,k) +! On choisit une relaxation quadratique. + gamma(ig,k)=gamma0(ig,k)*sqrt(dua(ig,k)**2+dva(ig,k)**2) + zua(ig,k)=(fm(ig,k)*zua(ig,k-1) & + & +(zf2*entr(ig,k)+gamma(ig,k))*zu(ig,k)) & + & /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)*zf*zf2 & + & +gamma(ig,k)) + zva(ig,k)=(fm(ig,k)*zva(ig,k-1) & + & +(zf2*entr(ig,k)+gamma(ig,k))*zv(ig,k)) & + & /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)*zf*zf2 & + & +gamma(ig,k)) + +! print*,' OUTPUT DV2 !!!!!!!!!!!!!!!!!!!!',k,zua(ig,k),zva(ig,k),zu(ig,k),zv(ig,k),dua(ig,k),dva(ig,k) + dua(ig,k)=zua(ig,k)-zu(ig,k) + dva(ig,k)=zva(ig,k)-zv(ig,k) + ue(ig,k)=(zu(ig,k)-zf*zua(ig,k))*zf2 + ve(ig,k)=(zv(ig,k)-zf*zva(ig,k))*zf2 + endif + enddo +! Fin des iterations +!-------------------- + enddo + + enddo ! k=2,nlay + +! Calcul du flux vertical de moment dans l'environnement. +!--------------------------------------------------------- + do k=2,nlay + do ig=1,ngrid + wud(ig,k)=fm(ig,k)*ue(ig,k) + wvd(ig,k)=fm(ig,k)*ve(ig,k) + enddo + enddo + do ig=1,ngrid + wud(ig,1)=0. + wud(ig,nlay+1)=0. + wvd(ig,1)=0. + wvd(ig,nlay+1)=0. + enddo + +! calcul des tendances. +!----------------------- + do k=1,nlay + do ig=1,ngrid + pduadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zua(ig,k) & + & -(entr(ig,k)+gamma(ig,k))*ue(ig,k) & + & -wud(ig,k)+wud(ig,k+1)) & + & /masse(ig,k) + pdvadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zva(ig,k) & + & -(entr(ig,k)+gamma(ig,k))*ve(ig,k) & + & -wvd(ig,k)+wvd(ig,k+1)) & + & /masse(ig,k) + enddo + enddo + + +! Sorties eventuelles. +!---------------------- + +! if (nlarga0>0) then +! print*,'WARNING !!!!!! DANS THERMCELL_DV2 ' +! print*,nlarga0,' points pour lesquels laraga=0. dans un thermique' +! print*,'Il faudrait decortiquer ces points' +! endif + +! =========================================================================== +! ============= FIN DV2 ===================================================== +! =========================================================================== + + else + + modname='momentum' + call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr, & + & masse,zu,pduadj,ztvd,fm_down,ztv,modname,lmax) + + call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr, & + & masse,zv,pdvadj,ztvd,fm_down,ztv,modname,lmax) + + endif + +!------------------------------------------------------------------ +! incrementation dt +!------------------------------------------------------------------ + + do l=1,nlay + do ig=1,ngrid + pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l) + enddo + enddo + +!------------------------------------------------------------------ +! calcul du transport vertical de traceurs +!------------------------------------------------------------------ + + if (nq .ne. 0.) then + modname='tracer' + DO iq=1,nq + call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr, & + & masse,pq(:,:,iq),pdqadj(:,:,iq),ztvd,fm_down,ztv,modname,lmax) + + ENDDO + endif + +!------------------------------------------------------------------ +! calcul du transport vertical de la tke +!------------------------------------------------------------------ + + modname='tke' + call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr, & + & masse,pq2,pdq2adj,ztvd,fm_down,ztv,modname,lmax) + +! =========================================================================== +! ============= FIN TRANSPORT =============================================== +! =========================================================================== + + +!------------------------------------------------------------------ +! Calculs de diagnostiques pour les sorties +!------------------------------------------------------------------ +! DIAGNOSTIQUE +! We compute interface values for teta env and th. The last interface +! value does not matter, as the mass flux is 0 there. + + + do l=1,nlay-1 + do ig=1,ngrid + teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l)) + teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l)) + teta_env_int(ig,l)=0.5*(ztv(ig,l+1)+ztv(ig,l)) + enddo + enddo + do ig=1,ngrid + teta_th_int(ig,nlay)=teta_th_int(ig,nlay-1) + teta_env_int(ig,nlay)=teta_env_int(ig,nlay-1) + teta_down_int(ig,nlay)=teta_down_int(ig,nlay-1) + enddo + do l=1,nlay + do ig=1,ngrid + heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l)) + buoyancyOut(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l) + buoyancyEst(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l) + heatFlux_down(ig,l)=fm_down(ig,l)*(teta_down_int(ig,l)-teta_env_int(ig,l))/rhobarz(ig,l) + enddo + enddo + + return + end Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 148) +++ trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 161) @@ -638,4 +638,13 @@ ENDIF + if (calltherm .and. outptherm) then + if (ngrid .eq. 1) then + call WRITEDIAGFI(ngrid,'zh','zh inside vdifc', + & 'SI',1,ph(:,:)+pdhfi(:,:)*ptimestep) + call WRITEDIAGFI(ngrid,'zkh','zkh', + & 'SI',1,zkh) + endif + endif + RETURN END