Index: trunk/LMDZ.MARS/libf/phymars/inifis.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/inifis.F (revision 565) +++ trunk/LMDZ.MARS/libf/phymars/inifis.F (revision 566) @@ -249,5 +249,5 @@ $ "1-> new correction", $ "(matters only if callnirco2=T)" - nircorr=0 !default value + nircorr=1 !default value call getin("nircorr",nircorr) write(*,*) " nircorr = ",nircorr @@ -290,5 +290,4 @@ endif -#ifndef MESOSCALE if (calladj .and. callrichsl .and. (.not. calltherm)) then print*,'You should not be calling the convective adjustment @@ -301,5 +300,5 @@ stop endif -#endif + write(*,*) "call CO2 condensation ?" callcond=.true. ! default value Index: trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F (revision 565) +++ trunk/LMDZ.MARS/libf/phymars/pbl_parameters.F (revision 566) @@ -1,4 +1,4 @@ SUBROUTINE pbl_parameters(ngrid,nlay,ps,pplay,pz0, - & pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,pts,ph,z_out, + & pg,zzlay,zzlev,pu,pv,wstar_in,hfmax,zmax,pts,ph,z_out,n_out, & Teta_out,u_out,ustar,tstar,L_mo,vhf,vvv) IMPLICIT NONE @@ -30,11 +30,12 @@ ! pts(ngrid) surface temperature ! ph(ngrid,nlay) potential temperature T*(p/ps)^kappa -! z_out height of interpolation +! z_out(n_out) heights of interpolation +! n_out number of points for interpolation ! ! outputs: ! ------ ! -! Teta_out(ngrid) interpolated teta -! u_out(ngrid) interpolated u +! Teta_out(ngrid,n_out) interpolated teta +! u_out(ngrid,n_out) interpolated u ! ustar(ngrid) friction velocity ! tstar(ngrid) friction temperature @@ -55,5 +56,5 @@ ! ---------- - INTEGER, INTENT(IN) :: ngrid,nlay + INTEGER, INTENT(IN) :: ngrid,nlay,n_out REAL, INTENT(IN) :: pz0(ngrid),ps(ngrid),pplay(ngrid,nlay) REAL, INTENT(IN) :: pg,zzlay(ngrid,nlay),zzlev(ngrid,nlay) @@ -61,11 +62,11 @@ REAL, INTENT(IN) :: wstar_in(ngrid),hfmax(ngrid),zmax(ngrid) REAL, INTENT(IN) :: pts(ngrid),ph(ngrid,nlay) - REAL, INTENT(IN) :: z_out + REAL, INTENT(IN) :: z_out(n_out) ! Outputs: ! -------- - REAL, INTENT(OUT) :: Teta_out(ngrid),u_out(ngrid) - REAL T_out(ngrid) + REAL, INTENT(OUT) :: Teta_out(ngrid,n_out),u_out(ngrid,n_out) + REAL T_out(ngrid,n_out) REAL, INTENT(OUT) :: ustar(ngrid), tstar(ngrid) REAL wstar(ngrid) @@ -75,5 +76,5 @@ ! ------ - INTEGER ig,k + INTEGER ig,k,n REAL karman,nu DATA karman,nu/.41,0.001/ @@ -118,4 +119,6 @@ !------------------------------------------------------------------------ + DO n=1,n_out + c Initialisation : @@ -123,5 +126,5 @@ ustar(:)=0. tstar(:)=0. - zout=z_out + zout=z_out(n) reynolds(:)=0. pz0t = 0. @@ -247,10 +250,10 @@ DO ig=1,ngrid IF(zout .lt. pz0tcomp(ig)) THEN - u_out(ig)=0. - Teta_out(ig)=pts(ig) + u_out(ig,n)=0. + Teta_out(ig,n)=pts(ig) ELSEIF (L_mo(ig) .gt. 0.) THEN - u_out(ig)=(ustar(ig)/karman)*log(zout/pz0(ig)) + + u_out(ig,n)=(ustar(ig)/karman)*log(zout/pz0(ig)) + & 5.*(ustar(ig)/(karman*L_mo(ig)))*(zout-pz0(ig)) - Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman)) + Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman)) & *log(zout/pz0tcomp(ig)) + & 5.*(tstar(ig)/(prandtl(ig)*karman*L_mo(ig))) @@ -260,5 +263,5 @@ IF(L_mo(ig) .gt. -1000.) THEN - u_out(ig)=(ustar(ig)/karman)*( + u_out(ig,n)=(ustar(ig)/karman)*( & 2.*atan((1.-16.*zout/L_mo(ig))**0.25) & -2.*atan((1.-16.*pz0(ig)/L_mo(ig))**0.25) @@ -270,5 +273,5 @@ & ) - Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( + Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( & 2.*log(1.+sqrt(1.-16.*pz0tcomp(ig)/L_mo(ig))) & -2.*log(1.+sqrt(1.-16.*zout/L_mo(ig))) @@ -285,5 +288,5 @@ ! (we do that to avoid using r*4 precision, otherwise, we get -inf values) - u_out(ig)=(ustar(ig)/karman)*( + u_out(ig,n)=(ustar(ig)/karman)*( & (4./L_mo(ig))*(zout-pz0(ig)) & + (20./(L_mo(ig))**2)*(zout**2-pz0(ig)**2) @@ -292,5 +295,5 @@ & ) - Teta_out(ig)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( + Teta_out(ig,n)=pts(ig)+(tstar(ig)/(prandtl(ig)*karman))*( & (8./L_mo(ig))*(zout-pz0tcomp(ig)) & + (48./(L_mo(ig))**2)*(zout**2-pz0tcomp(ig)**2) @@ -301,9 +304,9 @@ ENDIF ELSE - u_out(ig)=0. - Teta_out(ig)=pts(ig) + u_out(ig,n)=0. + Teta_out(ig,n)=pts(ig) ENDIF IF(zout .lt. pz0(ig)) THEN - u_out(ig)=0. + u_out(ig,n)=0. ENDIF ENDDO @@ -314,12 +317,13 @@ IF ((.not.calltherm).and.(calladj)) THEN - Teta_out(:)=ph(:,1) + Teta_out(:,n)=ph(:,1) ENDIF - T_out(:) = Teta_out(:)*(exp( + T_out(:,n) = Teta_out(:,n)*(exp( & (zout/zzlay(:,1))*(log(pplay(:,1)/ps)) & ) & )**rcp + ENDDO !of n=1,n_out !------------------------------------------------------------------------ @@ -333,4 +337,6 @@ ! Nearest index for the pbl height + + IF (calltherm) THEN pbl_height_index(:)=1 @@ -361,4 +367,5 @@ ENDDO +! Recompute wstar ! We follow Spiga et. al 2010 (QJRMS) ! ------------ @@ -410,46 +417,5 @@ vvv(:) = dvvv(:)*(wstar(:))**2 -!------------------------------------------------------------------------ -!------------------------------------------------------------------------ -! OUTPUTS -!------------------------------------------------------------------------ -!------------------------------------------------------------------------ - - IF (ngrid .eq. 1) THEN - dimout=0 - ELSE - dimout=2 - ENDIF - - call WRITEDIAGFI(ngrid,'Teta_out', - & 'potential temperature at z_out','K', - & dimout,Teta_out) - call WRITEDIAGFI(ngrid,'u_out', - & 'horizontal velocity norm at z_out','m/s', - & dimout,u_out) - call WRITEDIAGFI(ngrid,'u_star', - & 'friction velocity','m/s', - & dimout,ustar) - call WRITEDIAGFI(ngrid,'teta_star', - & 'friction potential temperature','K', - & dimout,tstar) - call WRITEDIAGFI(ngrid,'L', - & 'Monin Obukhov length','m', - & dimout,L_mo) -! call WRITEDIAGFI(ngrid,'w_star', -! & 'Free convection velocity','m', -! & dimout,wstar) -! call WRITEDIAGFI(ngrid,'z_i', -! & 'PBL height','m', -! & dimout,zmax) -! call WRITEDIAGFI(ngrid,'hf_max', -! & 'Maximum vertical heat flux','m', -! & dimout,hfmax) - call WRITEDIAGFI(ngrid,'vvv', - & 'Vertical velocity variance at zout','m', - & dimout,vvv) - call WRITEDIAGFI(ngrid,'vhf', - & 'Vertical heat flux at zout','m', - & dimout,vhf) + ENDIF ! of if calltherm RETURN Index: trunk/LMDZ.MARS/libf/phymars/physiq.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 565) +++ trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 566) @@ -345,9 +345,11 @@ REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx) REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx) - INTEGER lmax_th(ngridmx) + INTEGER lmax_th(ngridmx),dimout,n_out,n + CHARACTER(50) zstring REAL dtke_th(ngridmx,nlayermx+1) REAL zcdv(ngridmx), zcdh(ngridmx) - REAL Teta_out(ngridmx),u_out(ngridmx) ! Interpolated teta and u at z_out - REAL z_out ! height of interpolation between z0 and z1 [meters] + REAL, ALLOCATABLE, DIMENSION(:,:) :: Teta_out + REAL, ALLOCATABLE, DIMENSION(:,:) :: u_out ! Interpolated teta and u at z_out + REAL, ALLOCATABLE, DIMENSION(:) :: z_out ! height of interpolation between z0 and z1 [meters] REAL ustar(ngridmx),tstar(ngridmx) ! friction velocity and friction potential temp REAL L_mo(ngridmx),wstarpbl(ngridmx),vhf(ngridmx),vvv(ngridmx) @@ -838,8 +840,7 @@ c----------------------------------------------------------------------- -c TEST. Thermals : -c HIGHLY EXPERIMENTAL, BEWARE !! +c 5. Thermals : c ----------------------------- - + if(calltherm) then @@ -1377,4 +1378,60 @@ ENDIF ! of IF (lwrite) +c ---------------------------------------------------------- +c ---------------------------------------------------------- +c INTERPOLATIONS IN THE SURFACE-LAYER +c ---------------------------------------------------------- +c ---------------------------------------------------------- + + IF (1 .eq. 0.) THEN + IF (callrichsl) THEN + n_out=5 + + ALLOCATE(z_out(n_out)) + ALLOCATE(Teta_out(ngrid,n_out)) + ALLOCATE(u_out(ngrid,n_out)) + + z_out(:)=[0.001,0.05,0.1,0.5,1.] + u_out(:,:)=0. + Teta_out(:,:)=0. + + call pbl_parameters(ngrid,nlayer,ps,zplay,z0, + & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out,n_out, + & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) + +#ifndef MESOSCALE + IF (ngrid .eq. 1) THEN + dimout=0 + ELSE + dimout=2 + ENDIF + DO n=1,n_out + write(zstring, '(F9.6)') z_out(n) + call WRITEDIAGFI(ngrid,'Teta_out_'//trim(zstring), + & 'potential temperature at z_out','K',dimout,Teta_out(:,n)) + call WRITEDIAGFI(ngrid,'u_out_'//trim(zstring), + & 'horizontal velocity norm at z_out','m/s',dimout,u_out(:,n)) + ENDDO + call WRITEDIAGFI(ngrid,'u_star', + & 'friction velocity','m/s',dimout,ustar) + call WRITEDIAGFI(ngrid,'teta_star', + & 'friction potential temperature','K',dimout,tstar) + call WRITEDIAGFI(ngrid,'L', + & 'Monin Obukhov length','m',dimout,L_mo) + call WRITEDIAGFI(ngrid,'vvv', + & 'Vertical velocity variance at zout','m',dimout,vvv) + call WRITEDIAGFI(ngrid,'vhf', + & 'Vertical heat flux at zout','m',dimout,vhf) +#endif + + ENDIF + ENDIF ! of pbl interpolation outputs + +c ---------------------------------------------------------- +c ---------------------------------------------------------- +c END OF SURFACE LAYER INTERPOLATIONS +c ---------------------------------------------------------- +c ---------------------------------------------------------- + IF (ngrid.NE.1) THEN @@ -1863,25 +1920,4 @@ c ---------------------------------------------------------- - -c ---------------------------------------------------------- -c Outputs of surface layer -c ---------------------------------------------------------- - - - z_out=0. - if (calltherm .and. (z_out .gt. 0.)) then - - call pbl_parameters(ngrid,nlayer,ps,zplay,z0, - & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, - & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) - - else - if((.not. calltherm).and.(z_out .gt. 0.)) then - print*, 'WARNING : no interpolation in surface-layer :' - print*, 'Outputing surface-layer quantities without thermals - & does not make sense' - endif - endif - c ---------------------------------------------------------- c Outputs of thermals @@ -1961,20 +1997,4 @@ ! THERMALS STUFF 1D - - z_out=1. - if (calltherm .and. (z_out .gt. 0.)) then - - call pbl_parameters(ngrid,nlayer,ps,zplay,z0, - & g,zzlay,zzlev,zu,zv,wstar,hfmax_th,zmax_th,tsurf,zh,z_out, - & Teta_out,u_out,ustar,tstar,wstarpbl,L_mo,vhf,vvv) - - else - if((.not. calltherm).and.(z_out .gt. 0.)) then - print*, 'WARNING : no interpolation in surface-layer :' - print*, 'Outputing surface-layer quantities without thermals - & does not make sense' - endif - endif - if(calltherm) then