Index: /trunk/LMDZ.MARS/README =================================================================== --- /trunk/LMDZ.MARS/README (revision 336) +++ /trunk/LMDZ.MARS/README (revision 337) @@ -1120,2 +1120,42 @@ --> Added settings which works for 35 levels in inidissip.F according to FL. Commented for the moment. --> Checked compilation and run, looks fine. Note that 'jmars.20101220' is needed. + +== 10/10/2011 == AC + +*********** +This commit aims at increasing the thermals speed, especially for large tracer number configurations. The idea behind this commit is to advect non-active conserved variables outside of the sub-timestep of the thermals. Because these variables are not used in thermals computation, we can decouple them: + +momentum: can be decoupled because we assume a constant ratio between horizontal velocity in alimentation layer and maximum vertical velocity in the thermal +s. + +tracer: can be decoupled because we do not take condensation of any tracer into account and hence do not liberate latent heat nor form clouds in the thermal +s. + +temperature: cannot be decoupled (of course) +*********** + +D 336 libf/phymars/thermcell_dqupdown.F90 +^---------------- Deleted and replaced by a simpler version. Notes about downdraft advection are still available from revision 336 of SVN in thermcell_dqu +pdown. + +A 0 libf/phymars/thermcell_dqup.F90 +^---------------- New upward advection for tracers and momentum in thermals. Several changes are done compared to the old approach: + - Updraft quantities are not longer computed by making hypothesis on the amount of advected air. + - In general, the formalism for updraft computation is much simpler and clearer. + - Tracer tendancies are no longer computed using the conservation equation. Instead, we use the divergence + of an approximated turbulent flux of the concerned quantity, where downdraft are also neglected. + +M 336 libf/phymars/thermcell_main_mars.F90 +^---------------- The Main does not call anymore thermcell_dqupdown, which it was doing 2+tracer_number times per subtimestep (140 times per physical step for a 2 tracer config) + +M 336 libf/phymars/calltherm_mars.F90 +^---------------- Entrainment, detrainment and mass-fluxes are recomputed in the sub-timestep loop. Their final value after iterations is used by the new + advection routine to compute tracer and momentum fluxes. + +*********** Results + +- Conservation of tracers has been assessed over 1 yr in 1D and found to be comparable to that obtained with the simple convective adjustment. (it actually +seems to be better by a factor of 10%!) +- GCM speed-up is of about 20% compared to the old thermal configuration, for a 2 tracer case. +- Advection of sharp tracer profiles has been successfully observed, similar to the old method. + Index: /trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90 (revision 336) +++ /trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90 (revision 337) @@ -2,6 +2,6 @@ ! $Id: calltherm.F90 1428 2010-09-13 08:43:37Z fairhead $ ! - subroutine calltherm_mars(dtime,zzlev,zzlay & - & ,pplay,paprs,pphi & + subroutine calltherm_mars(ptimestep,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 & @@ -15,9 +15,10 @@ #include "dimensions.h" #include "dimphys.h" - - REAL dtime +#include "comcstfi.h" + + REAL ptimestep LOGICAL logexpr0, logexpr2(ngridmx,nlayermx), logexpr1(ngridmx) REAL fact - INTEGER nbptspb + INTEGER nbptspb,iq,l REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx) @@ -27,5 +28,5 @@ REAL t_seri(ngridmx,nlayermx),pq_therm(ngridmx,nlayermx,nqmx) REAL q2_therm(ngridmx,nlayermx) - REAL paprs(ngridmx,nlayermx+1) + REAL pplev(ngridmx,nlayermx+1) REAL pplay(ngridmx,nlayermx) REAL pphi(ngridmx,nlayermx) @@ -41,4 +42,5 @@ real fm_therm(ngridmx,nlayermx+1) real entr_therm(ngridmx,nlayermx),detr_therm(ngridmx,nlayermx) + REAL masse(ngridmx,nlayermx) !******************************************************** @@ -51,4 +53,6 @@ real lmax_real(ngridmx) real zmax(ngridmx),zmaxth(ngridmx) + REAL zdz(ngridmx,nlayermx) + !nouvelles variables pour la convection @@ -76,5 +80,5 @@ real zbuoyancyEst(ngridmx,nlayermx) - character (len=20) :: modname='calltherm' + character (len=20) :: modname character (len=80) :: abort_message @@ -107,7 +111,7 @@ call getin("r_aspect_thermals",r_aspect_thermals) -! fm_therm(:,:)=0. -! detr_therm(:,:)=0. -! entr_therm(:,:)=0. + fm_therm(:,:)=0. + detr_therm(:,:)=0. + entr_therm(:,:)=0. heatFlux(:,:)=0. @@ -120,5 +124,5 @@ lmax_real(:)=0. - zdt=dtime/REAL(nsplit_thermals) + zdt=ptimestep/REAL(nsplit_thermals) do isplit=1,nsplit_thermals @@ -130,7 +134,7 @@ ! cas de splitting -! zfm_therm(:,:)=0. -! zentr_therm(:,:)=0. -! zdetr_therm(:,:)=0. + zfm_therm(:,:)=0. + zentr_therm(:,:)=0. + zdetr_therm(:,:)=0. ! zheatFlux(:,:)=0. @@ -153,5 +157,5 @@ CALL thermcell_main_mars(zdt & ! CALL thermcell_main_mars_coupled_v2(zdt & - & ,pplay,paprs,pphi,zzlev,zzlay & + & ,pplay,pplev,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 & @@ -165,21 +169,21 @@ ! 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 + d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact +! d_u_the(:,:)=d_u_the(:,:)*fact +! d_v_the(:,:)=d_v_the(:,:)*fact ! dq2_the(:,:)=dq2_the(:,:)*fact - if (nqmx .ne. 0) then - d_q_the(:,:,:)=d_q_the(:,:,:)*fact - endif +! if (nqmx .ne. 0) then +! d_q_the(:,:,:)=d_q_the(:,:,:)*fact +! endif zmaxth(:)=zmaxth(:)+zmax(:)*fact lmax_real(:)=lmax_real(:)+float(lmax(:))*fact -! fm_therm(:,:)=fm_therm(:,:) & -! & +zfm_therm(:,:)*fact -! entr_therm(:,:)=entr_therm(:,:) & -! & +zentr_therm(:,:)*fact -! detr_therm(:,:)=detr_therm(:,:) & -! & +zdetr_therm(:,:)*fact + fm_therm(:,:)=fm_therm(:,:) & + & +zfm_therm(:,:)*fact + entr_therm(:,:)=entr_therm(:,:) & + & +zentr_therm(:,:)*fact + detr_therm(:,:)=detr_therm(:,:) & + & +zdetr_therm(:,:)*fact heatFlux(:,:)=heatFlux(:,:) & @@ -197,14 +201,14 @@ 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(:,:,:) +! 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(:,:,:) +! 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(:,:) @@ -218,17 +222,45 @@ !**************************************************************** -! do i=1,ngridmx -! do k=1,nlayermx -! if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0. -! print*,'youpi je sers a quelque chose !' -! enddo -! enddo - - DO i=1,ngridmx - hfmax(i)=MAXVAL(heatFlux(i,:)+heatFlux_down(i,:)) - wmax(i)=MAXVAL(zw2(i,:)) - ENDDO - - lmax(:)=nint(lmax_real(:)) +! Now that we have computed total entrainment and detrainment, we can +! advect u, v, and q in thermals. (theta already advected). We can do +! that separatly because u,v,and q are not used in thermcell_main for +! any thermals-related computation : they are purely passive. + +!calcul de la masse + do l=1,nlayermx + masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g + enddo + +!calcul de l'epaisseur des couches + do l=1,nlayermx + zdz(:,l)=zzlev(:,l+1)-zzlev(:,l) + enddo + + + modname='momentum' + call thermcell_dqup(ngridmx,nlayermx,ptimestep & + & ,fm_therm,entr_therm,detr_therm, & + & masse,u_seri,d_u_ajs,modname,zdz) + + call thermcell_dqup(ngridmx,nlayermx,ptimestep & + & ,fm_therm,entr_therm,detr_therm, & + & masse,v_seri,d_v_ajs,modname,zdz) + + if (nqmx .ne. 0.) then + modname='tracer' + DO iq=1,nqmx + call thermcell_dqup(ngridmx,nlayermx,ptimestep & + & ,fm_therm,entr_therm,detr_therm, & + & masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),modname,zdz) + + ENDDO + endif + + DO i=1,ngridmx + hfmax(i)=MAXVAL(heatFlux(i,:)+heatFlux_down(i,:)) + wmax(i)=MAXVAL(zw2(i,:)) + ENDDO + + lmax(:)=nint(lmax_real(:)) return Index: /trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90 (revision 337) +++ /trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90 (revision 337) @@ -0,0 +1,105 @@ + subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm0,entr0, & + & detr0,masse0,q_therm,dq_therm,charvar,zdz) + 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 +!======================================================================= + +#include "dimensions.h" +#include "dimphys.h" + +! ============================ INPUTS ============================ + + INTEGER, INTENT(IN) :: ngrid,nlayer + REAL, INTENT(IN) :: ptimestep + REAL, INTENT(IN) :: fm0(ngridmx,nlayermx+1) + REAL, INTENT(IN) ::entr0(ngridmx,nlayermx),detr0(ngridmx,nlayermx) + REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx) + CHARACTER (LEN=20), INTENT(IN) :: charvar + REAL, INTENT(IN) :: masse0(ngridmx,nlayermx) + REAL, INTENT(IN) :: zdz(ngridmx,nlayermx) + +! ============================ OUTPUTS =========================== + + REAL, INTENT(OUT) :: dq_therm(ngridmx,nlayermx) + +! ============================ LOCAL ============================= + + REAL q(ngridmx,nlayermx) + REAL qa(ngridmx,nlayermx) + REAL qd(ngridmx,nlayermx) + INTEGER ig,k + REAL gammaf(ngridmx,nlayermx) + +! =========== Init ============================================== + + qa(:,:)=q_therm(:,:) + q(:,:)=q_therm(:,:) + + do ig=1,ngridmx + do k=1, nlayermx + if (fm0(ig,k)+entr0(ig,k) .gt. 0.) then + gammaf(ig,k)=fm0(ig,k)/(fm0(ig,k)+entr0(ig,k)) + else + gammaf(ig,k)=0. + endif + enddo + enddo + + +! =========== Updraft ============================================ + +! qa(:,1)=q(:,1) + + do ig=1,ngridmx + do k=2, nlayermx + + qa(ig,k)=gammaf(ig,k)*qa(ig,k-1) +(1.-gammaf(ig,k))*q(ig,k) + + enddo + enddo + + + +! ====== dq ====================================================== + +! do ig=1,ngridmx +! 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) +! enddo +! do k=2,nlayermx-1 +! do ig=1, ngridmx +! 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) +! enddo +! enddo +! +! do ig=1, ngridmx +! dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx) & +! & -entr0(ig,nlayermx)*q(ig,nlayermx) & +! & -fm0(ig,nlayermx)*q(ig,nlayermx)) & +! & *ptimestep/masse0(ig,nlayermx) +! +! enddo + + do ig=1, ngridmx + do k=1,nlayermx-1 + dq_therm(ig,k)=-(ptimestep/masse0(ig,k))*( & + & fm0(ig,k+1)*(qa(ig,k+1)-q(ig,k+1)) - & + & fm0(ig,k)*(qa(ig,k)-q(ig,k)) ) & + & /zdz(ig,k) + + enddo + enddo + + return + end Index: /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90 (revision 336) +++ /trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90 (revision 337) @@ -38,7 +38,7 @@ REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx) - REAL, INTENT(OUT) :: pduadj(ngridmx,nlayermx) - REAL, INTENT(OUT) :: pdvadj(ngridmx,nlayermx) - REAL, INTENT(OUT) :: pdqadj(ngridmx,nlayermx,nqmx) + REAL :: pduadj(ngridmx,nlayermx) + REAL :: pdvadj(ngridmx,nlayermx) + REAL :: pdqadj(ngridmx,nlayermx,nqmx) ! REAL, INTENT(OUT) :: pdq2adj(ngridmx,nlayermx) REAL :: pdq2adj(ngridmx,nlayermx) @@ -185,6 +185,6 @@ detr(:,:)=0. fm(:,:)=0. - zu(:,:)=pu(:,:) - zv(:,:)=pv(:,:) +! zu(:,:)=pu(:,:) +! zv(:,:)=pv(:,:) ztv(:,:)=pt(:,:)/zpopsk(:,:) @@ -1308,5 +1308,5 @@ ! 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) + 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)) & @@ -1377,10 +1377,10 @@ else - modname='momentum' - call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, & - & masse,zu,pduadj,ztvd,fm_down,ztv,modname,lmax) - - call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, & - & masse,zv,pdvadj,ztvd,fm_down,ztv,modname,lmax) +! modname='momentum' +! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, & +! & masse,zu,pduadj,ztvd,fm_down,ztv,modname,lmax) +! +! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, & +! & masse,zv,pdvadj,ztvd,fm_down,ztv,modname,lmax) endif @@ -1400,12 +1400,12 @@ !------------------------------------------------------------------ - if (nqmx .ne. 0.) then - modname='tracer' - DO iq=1,nqmx - call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, & - & masse,pq(:,:,iq),pdqadj(:,:,iq),ztvd,fm_down,ztv,modname,lmax) - - ENDDO - endif +! if (nqmx .ne. 0.) then +! modname='tracer' +! DO iq=1,nqmx +! call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr, & +! & masse,pq(:,:,iq),pdqadj(:,:,iq),ztvd,fm_down,ztv,modname,lmax) +! +! ENDDO +! endif !------------------------------------------------------------------