Index: trunk/LMDZ.MARS/libf/phymars/callkeys.h =================================================================== --- trunk/LMDZ.MARS/libf/phymars/callkeys.h (revision 283) +++ trunk/LMDZ.MARS/libf/phymars/callkeys.h (revision 284) @@ -12,5 +12,5 @@ & ,callg2d,linear,rayleigh,tracer,active,doubleq,submicron & & ,lifting,callddevil,scavenging,sedimentation,activice,water & - & ,caps,photochem,calltherm,outptherm,callslope + & ,caps,photochem,calltherm,outptherm,callrichsl,callslope COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche & @@ -24,5 +24,5 @@ & ,callnirco2,callnlte,callthermos,callconduct, & & calleuv,callmolvis,callmoldiff,thermochem,thermoswater & - & ,calltherm,outptherm,callslope + & ,calltherm,outptherm,callrichsl,callslope Index: trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90 (revision 283) +++ trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90 (revision 284) @@ -71,6 +71,4 @@ REAL buoyancyEst(ngridmx,nlayermx) REAL hfmax(ngridmx),wmax(ngridmx) - - REAL tstart,tstop !--------------------------------------------------------- @@ -135,6 +133,4 @@ endif - call cpu_time(tstart) - CALL calltherm_mars(ptimestep,zzlev,zzlay & & ,pplay,pplev,pphi & @@ -146,7 +142,4 @@ & ,zpopsk,ztla,heatFlux,heatFlux_down & & ,buoyancyOut,buoyancyEst,hfmax,wmax) - call cpu_time(tstop) - print*,'TOTAL elapsed time in thermals : ',tstop-tstart - ! Accumulation des tendances. On n'accumule pas les quantités de traceurs car celle ci n'a pas du changer Index: trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90 =================================================================== --- trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90 (revision 283) +++ trunk/LMDZ.MARS/libf/phymars/calltherm_mars.F90 (revision 284) @@ -99,5 +99,5 @@ r_aspect_thermals=0.7 #else - nsplit_thermals=40 + nsplit_thermals=50 r_aspect_thermals=2. #endif Index: trunk/LMDZ.MARS/libf/phymars/inifis.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/inifis.F (revision 283) +++ trunk/LMDZ.MARS/libf/phymars/inifis.F (revision 284) @@ -252,4 +252,9 @@ stop endif + + write(*,*) "call Richardson-based surface layer ?" + callrichsl=.false. ! default value + call getin("callrichsl",callrichsl) + write(*,*) " callrichsl = ",callrichsl write(*,*) "call CO2 condensation ?" Index: trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F (revision 283) +++ trunk/LMDZ.MARS/libf/phymars/meso_inc/meso_inc_les.F (revision 284) @@ -1,21 +1,36 @@ - DO ig=1,ngrid - !! sensible heat flux in W/m2 -! hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig) + if (callrichsl .eq. .false.) then + + DO ig=1,ngrid +!! sensible heat flux in W/m2 + + hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig) + +!! u star in similarity theory in m/s + ust(ig) = 0.4 + . * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ) + . / log( 1.E+0 + zzlay(ig,1)/z0_default ) + ENDDO + + else + + DO ig=1,ngrid ! New SL parametrization, correct formulation for hfx : - hfx(ig) = (pplay(ig,1)/(r*pt(ig,1)))*cpp - & *sqrt((pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1))) - & *zcdh(ig)*(tsurf(ig)-zh(ig,1)) + hfx(ig) = (pplay(ig,1)/(r*pt(ig,1)))*cpp + & *sqrt(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + & + (1.0*wmax_th(ig))**2) + & *zcdh(ig)*(tsurf(ig)-zh(ig,1)) - !! u star in similarity theory in m/s -! ust(ig) = 0.4 -! . * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ) -! . / log( 1.E+0 + zzlay(ig,1)/z0_default ) ! New SL parametrization, ust is more accurately computed in vdif_cd : - ust(ig) = sqrt(zcdv(ig)*(pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1))) + ust(ig) = sqrt(zcdv(ig)* + & (pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) + (1.0*wmax_th(ig))**2) + & ) - ENDDO + ENDDO + + endif !of if callrichsl + ! write (*,*) 'PHYS HFX cp zdts', hfx(100), zflubid(100), ! . capcal(100), Index: trunk/LMDZ.MARS/libf/phymars/physiq.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 283) +++ trunk/LMDZ.MARS/libf/phymars/physiq.F (revision 284) @@ -694,4 +694,30 @@ enddo +c Treatment of a special case : using new surface layer (Richardson based) +c without using the thermals in gcm and mesoscale can yield problems in +c weakly unstable situations when winds are near to 0. For those cases, we add +c a unit subgrid gustiness. Remember that thermals should be used we using the +c Richardson based surface layer model. + +#ifdef MESOSCALE + IF (flag_LES .eq. .false.) THEN + IF ((calltherm .eq. .false.) .and. (callrichsl .eq. .true.)) THEN + DO ig=1, ngridmx + IF (zh(ig,1) .lt. tsurf(ig)) THEN + wmax_th(ig)=1. + ENDIF + ENDDO + ENDIF + ENDIF +#else + IF ((calltherm .eq. .false.) .and. (callrichsl .eq. .true.)) THEN + DO ig=1, ngridmx + IF (zh(ig,1) .lt. tsurf(ig)) THEN + wmax_th(ig)=1. + ENDIF + ENDDO + ENDIF +#endif + c Calling vdif (Martian version WITH CO2 condensation) CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk, @@ -1433,14 +1459,15 @@ call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2, & co2ice) + c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7", c & "K",2,zt(1,7)) - call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2, - & fluxsurf_lw) - call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2, - & fluxsurf_sw_tot) - call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2, - & fluxtop_lw) - call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2, - & fluxtop_sw_tot) +c call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2, +c & fluxsurf_lw) +c call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2, +c & fluxsurf_sw_tot) +c call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2, +c & fluxtop_lw) +c call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2, +c & fluxtop_sw_tot) call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) @@ -1734,8 +1761,4 @@ endif -! --- - - - if(calltherm) then Index: trunk/LMDZ.MARS/libf/phymars/vdif_cd.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/vdif_cd.F (revision 283) +++ trunk/LMDZ.MARS/libf/phymars/vdif_cd.F (revision 284) @@ -36,4 +36,5 @@ #include "comcstfi.h" +#include "callkeys.h" c Arguments: @@ -87,17 +88,20 @@ c ------------------ - reynolds(:)=0. - c Original formulation : -c DO ig=1,ngrid -c z1=1.E+0 + pz(ig,1)/pz0(ig) -c zcd0=karman/log(z1) -c zcd0=zcd0*zcd0 -c pcdv(ig)=zcd0 -c pcdh(ig)=zcd0 -c ENDDO + if(callrichsl .eq. .false.) then + + DO ig=1,ngrid + z1=1.E+0 + pz(ig,1)/pz0(ig) + zcd0=karman/log(z1) + zcd0=zcd0*zcd0 + pcdv(ig)=zcd0 + pcdh(ig)=zcd0 + ENDDO ! print*,'old : cd,ch; ',pcdv,pcdh + else + + reynolds(:)=0. c New formulation (AC) : @@ -224,5 +228,5 @@ ! Some useful diagnostics : -! call WRITEDIAGFI(ngrid,'RiB', +! call WRITEDIAGFI(ngrid,'RiB', ! & 'Bulk Richardson nb','no units', ! & 2,rib) @@ -241,4 +245,6 @@ + endif !of if call richardson surface layer + RETURN END Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 283) +++ trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 284) @@ -86,5 +86,5 @@ REAL zc(ngridmx,nlayermx),zd(ngridmx,nlayermx) REAL zcst1 - REAL zu2 + REAL zu2(ngridmx) EXTERNAL SSUM,SCOPY @@ -234,14 +234,16 @@ & ,zcdv_true,zcdh_true) - DO ig=1,ngrid - - zu2=pu(ig,1)*pu(ig,1)+pv(ig,1)*pv(ig,1) - & +(1.*wmax(ig))**2 ! subgrid gustiness added by quadratic interpolation with a coeff beta, here beta=1. (tuned from + + IF (callrichsl .eq. .true.) then + zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1) + & +(1.*wmax(:))**2 ! subgrid gustiness added by quadratic interpolation with a coeff beta, here beta=1. (tuned from ! LES comparisons. This parameter is linked to the thermals settings) - - zcdv(ig)=zcdv_true(ig)*sqrt(zu2) ! 1 / bulk aerodynamic momentum conductance - zcdh(ig)=zcdh_true(ig)*sqrt(zu2) ! 1 / bulk aerodynamic heat conductance + ELSE + zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1) + ENDIF + + zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance + zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance ! these are the quantities to be looked at when comparing surface layers of different models - ENDDO ! Some usefull diagnostics for the new surface layer parametrization : @@ -254,6 +256,6 @@ ! & 2,zcdh_true) ! call WRITEDIAGFI(ngridmx,'u*', -! & 'friction velocity','m/s', -! & 2,sqrt(zcdv_true(:))*sqrt(zu2)) +! & 'friction velocity','m/s', +! & 2,sqrt(zcdv_true(:))*sqrt(zu2(:))) ! call WRITEDIAGFI(ngridmx,'T*', ! & 'friction temperature','K', @@ -479,7 +481,7 @@ c ---------------------------------------------------------------- DO ig=1,ngrid - zu2=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1) - zcdv(ig)=zcdv_true(ig)*sqrt(zu2) - zcdh(ig)=zcdh_true(ig)*sqrt(zu2) + zu2(ig)=zu(ig,1)*zu(ig,1)+zv(ig,1)*zv(ig,1) + zcdv(ig)=zcdv_true(ig)*sqrt(zu2(ig)) + zcdh(ig)=zcdh_true(ig)*sqrt(zu2(ig)) ENDDO