Index: LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90 (revision 4653) +++ LMDZ6/trunk/libf/phylmd/coef_diff_turb_mod.F90 (revision 4654) @@ -167,8 +167,10 @@ ! iflag_pbl peut etre utilise comme longuer de melange IF (iflag_pbl.GE.31) THEN - CALL vdif_kcay(knon,dtime,RG,RD,ypaprs,yt, & + if (knon>0) then + CALL vdif_kcay(knon,klev,dtime,RG,RD,ypaprs,yt, & yzlev,yzlay,yu,yv,yteta, & ycdragm,yq2,q2diag,ykmm,ykmn,yustar, & iflag_pbl) + endif ELSE IF (iflag_pbl<20) THEN CALL yamada4(ni,nsrf,knon,dtime,RG,RD,ypaprs,yt, & Index: LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90 (revision 4653) +++ LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90 (revision 4654) @@ -1,5 +1,5 @@ -module lscp_ini_mod - -implicit none +MODULE lscp_ini_mod + +IMPLICIT NONE ! PARAMETERS for lscp: @@ -259,5 +259,7 @@ CALL cloudth_ini(iflag_cloudth_vert,iflag_ratqs) -end subroutine lscp_ini - -end module lscp_ini_mod +RETURN + +END SUBROUTINE lscp_ini + +END MODULE lscp_ini_mod Index: LMDZ6/trunk/libf/phylmd/lscp_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lscp_mod.F90 (revision 4653) +++ LMDZ6/trunk/libf/phylmd/lscp_mod.F90 (revision 4654) @@ -7,5 +7,5 @@ !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE lscp(klon,klev,dtime,missing_val, & - paprs,pplay,t,q,ptconv,ratqs, & + paprs,pplay,temp,q,ptconv,ratqs, & d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri, & pfraclr,pfracld, & @@ -91,4 +91,6 @@ ! USE print_control_mod, ONLY: prt_level, lunout + +! USE de modules contenant des fonctions. USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc USE lscp_tools_mod, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat @@ -96,4 +98,5 @@ USE ice_sursat_mod, ONLY : ice_sursat +! Use du module lscp_ini_mod contenant les constantes USE lscp_ini_mod, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc @@ -119,5 +122,5 @@ REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa] REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa] - REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K) + REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature (K) REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! total specific humidity (= vapor phase) [kg/kg] INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds @@ -194,68 +197,75 @@ !---------------- - REAL qsl(klon), qsi(klon) - REAL zct, zcl,zexpo - REAL ctot(klon,klev) - REAL ctot_vol(klon,klev) - REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 - REAL zdqsdT_raw(klon) - REAL gammasat(klon),dgammasatdt(klon) ! coefficient to make cold condensation at the correct RH and derivative wrt T - REAL Tbef(klon),qlbef(klon),DT(klon),num,denom - REAL cste - REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) - REAL Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) - REAL erf - REAL qcloud(klon), icefrac_mpc(klon,klev), qincloud_mpc(klon) - REAL zrfl(klon), zrfln(klon), zqev, zqevt - REAL zifl(klon), zifln(klon), ziflprev(klon),zqev0,zqevi, zqevti - REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon) - REAL zoliql(klon), zoliqi(klon) - REAL zt(klon) - REAL zdz(klon),zrho(klon,klev),iwc(klon) - REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon) - REAL zmelt,zrain,zsnow,zprecip - REAL dzfice(klon) - REAL zsolid - REAL qtot(klon), qzero(klon) - REAL dqsl(klon),dqsi(klon) - REAL smallestreal + REAL,DIMENSION(klon) :: qsl, qsi + REAL :: zct, zcl,zexpo + REAL, DIMENSION(klon,klev) :: ctot + REAL, DIMENSION(klon,klev) :: ctot_vol + REAL, DIMENSION(klon) :: zqs, zdqs + REAL :: zdelta, zcor, zcvm5 + REAL, DIMENSION(klon) :: zdqsdT_raw + REAL, DIMENSION(klon) :: gammasat,dgammasatdt ! coefficient to make cold condensation at the correct RH and derivative wrt T + REAL, DIMENSION(klon) :: Tbef,qlbef,DT + REAL :: num,denom + REAL :: cste + REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta + REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2 + REAL :: erf + REAL, DIMENSION(klon,klev) :: icefrac_mpc + REAL, DIMENSION(klon) :: qcloud, qincloud_mpc + REAL, DIMENSION(klon) :: zrfl, zrfln + REAL :: zqev, zqevt + REAL, DIMENSION(klon) :: zifl, zifln, ziflprev + REAL :: zqev0,zqevi, zqevti + REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn + REAL, DIMENSION(klon) :: zoliql, zoliqi + REAL, DIMENSION(klon) :: zt + REAL, DIMENSION(klon,klev) :: zrho + REAL, DIMENSION(klon) :: zdz,iwc + REAL :: zchau,zfroi + REAL, DIMENSION(klon) :: zfice,zneb,znebprecip + REAL :: zmelt,zrain,zsnow,zprecip + REAL, DIMENSION(klon) :: dzfice + REAL :: zsolid + REAL, DIMENSION(klon) :: qtot, qzero + REAL, DIMENSION(klon) :: dqsl,dqsi + REAL :: smallestreal ! Variables for Bergeron process - REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl - REAL zqpreci(klon), zqprecl(klon) + REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl + REAL, DIMENSION(klon) :: zqpreci, zqprecl ! Variables precipitation energy conservation - REAL zmqc(klon) - REAL zalpha_tr - REAL zfrac_lessi - REAL zprec_cond(klon) - REAL zmair(klon), zcpair, zcpeau - REAL zlh_solid(klon), zm_solid ! for liquid -> solid conversion - REAL zrflclr(klon), zrflcld(klon) - REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon) - REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon) - REAL ziflclr(klon), ziflcld(klon) - REAL znebprecipclr(klon), znebprecipcld(klon) - REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon) - REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon) - REAL velo(klon,klev), vr, ffallv - REAL qlmpc, qimpc, rnebmpc - REAL radocondi(klon,klev), radocondl(klon,klev) - REAL effective_zneb - REAL distcltop1D(klon), temp_cltop1D(klon) - - - INTEGER i, k, n, kk - INTEGER n_i(klon), iter + REAL, DIMENSION(klon) :: zmqc + REAL :: zalpha_tr + REAL :: zfrac_lessi + REAL, DIMENSION(klon) :: zprec_cond + REAL, DIMENSION(klon) :: zmair + REAL :: zcpair, zcpeau + REAL, DIMENSION(klon) :: zlh_solid + REAL :: zm_solid ! for liquid -> solid conversion + REAL, DIMENSION(klon) :: zrflclr, zrflcld + REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld + REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr + REAL, DIMENSION(klon) :: ziflclr, ziflcld + REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld + REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb + REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr + REAL, DIMENSION(klon,klev) :: velo + REAL :: vr, ffallv + REAL :: qlmpc, qimpc, rnebmpc + REAL, DIMENSION(klon,klev) :: radocondi, radocondl + REAL :: effective_zneb + REAL, DIMENSION(klon) :: distcltop1D, temp_cltop1D + + + INTEGER i, k, n, kk, iter + INTEGER, DIMENSION(klon) :: n_i INTEGER ncoreczq - INTEGER mpc_bl_points(klon,klev) - - LOGICAL lognormale(klon) - LOGICAL keepgoing(klon) - LOGICAL,SAVE :: appel1er - !$OMP THREADPRIVATE(appel1er) + INTEGER, DIMENSION(klon,klev) :: mpc_bl_points + + LOGICAL, DIMENSION(klon) :: lognormale + LOGICAL, DIMENSION(klon) :: keepgoing CHARACTER (len = 20) :: modname = 'lscp' CHARACTER (len = 80) :: abort_message - DATA appel1er /.TRUE./ !=============================================================================== @@ -264,4 +274,5 @@ ! Few initial checks + IF (iflag_fisrtilp_qsat .LT. 0) THEN @@ -332,4 +343,6 @@ fcontrN(:,:) = 0.0 fcontrP(:,:) = 0.0 +distcltop(:,:)=0. +temp_cltop(:,:)=0. !c_iso: variable initialisation for iso @@ -346,5 +359,5 @@ ! Initialisation temperature and specific humidity DO i = 1, klon - zt(i)=t(i,k) + zt(i)=temp(i,k) zq(i)=q(i,k) !c_iso init of iso @@ -382,5 +395,5 @@ zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i) ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer - zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) & + zt(i) = ( (temp(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) & / (zcpair + zmqc(i)*zcpeau) @@ -625,5 +638,5 @@ zq,zqta,fraca, & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & - ratqs,zqs,t, & + ratqs,zqs,temp, & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) @@ -637,5 +650,5 @@ zq,zqta,fraca, & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & - ratqs,zqs,t, & + ratqs,zqs,temp, & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) @@ -648,5 +661,5 @@ zq,zqta,fraca, & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & - ratqs,zqs,t, & + ratqs,zqs,temp, & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) @@ -658,5 +671,5 @@ CALL cloudth_mpc(klon,klev,k,iflag_mpc_bl,mpc_bl_points, & - t,ztv,zq,zqta,fraca, zpspsk, & + temp,ztv,zq,zqta,fraca, zpspsk, & paprs,pplay,ztla,zthl,ratqs,zqs,psfl, & qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol, & @@ -739,5 +752,5 @@ IF (iflag_t_glace.GE.4) THEN ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top - CALL distance_to_cloud_top(klon,klev,k,t,pplay,paprs,rneb,distcltop1D,temp_cltop1D) + CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D) ENDIF CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, distcltop1D(:),temp_cltop1D(:),zfice(:),dzfice(:)) @@ -780,5 +793,5 @@ !------------------------------------ - CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), & + CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, temp(i,k), zq(i), & gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), & rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), & @@ -837,5 +850,5 @@ ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top IF (iflag_t_glace.GE.4) THEN - CALL distance_to_cloud_top(klon,klev,k,t,pplay,paprs,rneb,distcltop1D,temp_cltop1D) + CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D) distcltop(:,k)=distcltop1D(:) temp_cltop(:,k)=temp_cltop1D(:) @@ -1196,5 +1209,5 @@ d_q(i,k) = zq(i) - q(i,k) ! c_iso: same for isotopes - d_t(i,k) = zt(i) - t(i,k) + d_t(i,k) = zt(i) - temp(i,k) ENDDO @@ -1245,5 +1258,5 @@ IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN - IF (t(i,k) .GE. t_glace_min) THEN + IF (temp(i,k) .GE. t_glace_min) THEN zalpha_tr = a_tr_sca(3) ELSE @@ -1269,5 +1282,5 @@ IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN - IF (t(i,kk) .GE. t_glace_min) THEN + IF (temp(i,kk) .GE. t_glace_min) THEN zalpha_tr = a_tr_sca(1) ELSE @@ -1318,4 +1331,7 @@ ENDIF + +RETURN + END SUBROUTINE lscp !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Index: LMDZ6/trunk/libf/phylmd/physiqex_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiqex_mod.F90 (revision 4653) +++ LMDZ6/trunk/libf/phylmd/physiqex_mod.F90 (revision 4654) @@ -9,18 +9,32 @@ & debut,lafin,pdtphys, & & paprs,pplay,pphi,pphis,presnivs, & - & u,v,rot,t,qx, & + & u,v,rot,temp,qx, & & flxmass_w, & & d_u, d_v, d_t, d_qx, d_ps) - USE dimphy, only : klon,klev - USE infotrac_phy, only : nqtot - USE geometry_mod, only : latitude -! USE comcstphy, only : rg - USE ioipsl, only : ymds2ju - USE phys_state_var_mod, only : phys_state_var_init - USE phyetat0_mod, only: phyetat0 - USE output_physiqex_mod, ONLY: output_physiqex + +USE dimphy, only : klon,klev +USE infotrac_phy, only : nqtot +USE geometry_mod, only : latitude +USE ioipsl, only : ymds2ju +USE phys_state_var_mod, only : phys_state_var_init +USE phyetat0_mod, only: phyetat0 +USE output_physiqex_mod, ONLY: output_physiqex +use vdif_ini, only : vdif_ini_ +USE lmdz_thermcell_ini, ONLY : thermcell_ini +USE ioipsl_getin_p_mod, ONLY : getin_p +USE wxios, ONLY: missing_val, using_xios +USE lscp_mod, ONLY : lscp +USE lscp_ini_mod, ONLY : lscp_ini +USE add_phys_tend_mod, ONLY : add_phys_tend + IMPLICIT none + +include "YOETHF.h" + + + + ! ! Routine argument: @@ -32,20 +46,110 @@ logical,intent(in) :: lafin ! signals last call to physics real,intent(in) :: pdtphys ! physics time step (s) - real,intent(in) :: paprs(klon,klev+1) ! interlayer pressure (Pa) - real,intent(in) :: pplay(klon,klev) ! mid-layer pressure (Pa) - real,intent(in) :: pphi(klon,klev) ! geopotential at mid-layer - real,intent(in) :: pphis(klon) ! surface geopotential - real,intent(in) :: presnivs(klev) ! pseudo-pressure (Pa) of mid-layers - real,intent(in) :: u(klon,klev) ! eastward zonal wind (m/s) - real,intent(in) :: v(klon,klev) ! northward meridional wind (m/s) - real,intent(in) :: rot(klon,klev) ! northward meridional wind (m/s) - real,intent(in) :: t(klon,klev) ! temperature (K) - real,intent(in) :: qx(klon,klev,nqtot) ! tracers (.../kg_air) - real,intent(in) :: flxmass_w(klon,klev) ! vertical mass flux - real,intent(out) :: d_u(klon,klev) ! physics tendency on u (m/s/s) - real,intent(out) :: d_v(klon,klev) ! physics tendency on v (m/s/s) - real,intent(out) :: d_t(klon,klev) ! physics tendency on t (K/s) - real,intent(out) :: d_qx(klon,klev,nqtot) ! physics tendency on tracers - real,intent(out) :: d_ps(klon) ! physics tendency on surface pressure + real,dimension(klon,klev+1),intent(in) :: paprs ! interlayer pressure (Pa) + real,dimension(klon,klev),intent(in) :: pplay ! mid-layer pressure (Pa) + real,dimension(klon,klev),intent(in) :: pphi ! geopotential at mid-layer + real,dimension(klon),intent(in) :: pphis ! surface geopotential + real,dimension(klev),intent(in) :: presnivs ! pseudo-pressure (Pa) of mid-layers + real,dimension(klon,klev),intent(in) :: u ! eastward zonal wind (m/s) + real,dimension(klon,klev),intent(in) :: v ! northward meridional wind (m/s) + real,dimension(klon,klev),intent(in) :: rot ! northward meridional wind (m/s) + real,dimension(klon,klev),intent(in) :: temp ! temperature (K) + real,dimension(klon,klev,nqtot),intent(in) :: qx ! tracers (.../kg_air) + real,dimension(klon,klev),intent(in) :: flxmass_w ! vertical mass flux + real,dimension(klon,klev),intent(out) :: d_u ! physics tendency on u (m/s/s) + real,dimension(klon,klev),intent(out) :: d_v ! physics tendency on v (m/s/s) + real,dimension(klon,klev),intent(out) :: d_t ! physics tendency on t (K/s) + real,dimension(klon,klev,nqtot),intent(out) :: d_qx ! physics tendency on tracers + real,dimension(klon),intent(out) :: d_ps ! physics tendency on surface pressure + + real, dimension(klon,klev) :: u_loc + real, dimension(klon,klev) :: v_loc + real, dimension(klon,klev) :: t_loc + real, dimension(klon,klev) :: h_loc + real, dimension(klon,klev) :: d_u_loc,d_v_loc,d_t_loc,d_h_loc + + real, dimension(klon,klev) :: d_u_dyn,d_v_dyn,d_t_dyn + real, dimension(klon,klev,nqtot) :: d_q_dyn + real, allocatable, dimension(:,:), save :: u_prev,v_prev,t_prev + real, allocatable, dimension(:,:,:), save :: q_prev +!$OMP THREADPRIVATE(u_prev,v_prev,t_prev,q_prev) + + + + real, dimension(klon,klev) :: d_u_vdif,d_v_vdif,d_t_vdif,d_h_vdif + real, dimension(klon,klev) :: d_u_the,d_v_the,d_t_the + real, dimension(klon,klev,nqtot) :: q_loc,d_q_loc,d_q_vdif,d_q_the + +real, dimension(klon) :: capcal,z0m,z0h,dtsrf,emis,fluxsrf,cdh,cdv,tsrf_ +real, dimension(klon,klev) :: zzlay,masse,zpopsk +real, dimension(klon,klev+1) :: zzlev,kz_v,kz_h,richardson + +real, save, allocatable, dimension(:) :: tsrf,f0,zmax0 +real, save, allocatable, dimension(:,:) :: q2 +!$OMP THREADPRIVATE(tsrf,q2,f0,zmax0) + + real,save :: ratqsbas=0.002,ratqshaut=0.3,ratqsp0=50000.,ratqsdp=20000. + !$OMP THREADPRIVATE(ratqsbas,ratqshaut,ratqsp0,ratqsdp) + + +real :: z1,z2,tau_thermals +logical :: lwrite +integer :: iflag_replay + +integer :: iflag_thermals=18 + +!-------------------------------------------------------------- +! Declaration lscp +!-------------------------------------------------------------- + INTEGER :: iflag_cld_th ! flag that determines the distribution of convective clouds ! IN + INTEGER :: iflag_ice_thermo! flag to activate the ice thermodynamics ! IN + LOGICAL :: ok_ice_sursat ! flag to determine if ice sursaturation is activated ! IN + LOGICAL, DIMENSION(klon,klev) :: ptconv ! grid points where deep convection scheme is active ! IN + REAL, DIMENSION(klon,klev) :: ztv ! virtual potential temperature [K] ! IN + REAL, DIMENSION(klon,klev) :: zqta ! specific humidity within thermals [kg/kg] ! IN + REAL, DIMENSION(klon,klev+1) :: frac_the,fm_the + REAL, DIMENSION(klon,klev) :: zpspsk ! exner potential (p/100000)**(R/cp) ! IN + REAL, DIMENSION(klon,klev) :: ztla ! liquid temperature within thermals [K] ! IN + REAL, DIMENSION(klon,klev) :: zthl ! liquid potential temperature [K] ! INOUT + REAL, DIMENSION(klon,klev) :: ratqs ! function of pressure that sets the large-scale ! INOUT + REAL, DIMENSION(klon,klev) :: beta ! conversion rate of condensed water ! INOUT + REAL, DIMENSION(klon,klev) :: rneb_seri ! fraction nuageuse en memoire ! INOUT + REAL, DIMENSION(klon,klev) :: d_t_lscp ! temperature increment [K] ! OUT + REAL, DIMENSION(klon,klev) :: d_q_lscp ! specific humidity increment [kg/kg] ! OUT + REAL, DIMENSION(klon,klev) :: d_ql_lscp ! liquid water increment [kg/kg] ! OUT + REAL, DIMENSION(klon,klev) :: d_qi_lscp ! cloud ice mass increment [kg/kg] ! OUT + REAL, DIMENSION(klon,klev) :: rneb ! cloud fraction [-] ! OUT + REAL, DIMENSION(klon,klev) :: rneblsvol ! cloud fraction per unit volume [-] ! OUT + REAL, DIMENSION(klon,klev) :: pfraclr ! precip fraction clear-sky part [-] ! OUT + REAL, DIMENSION(klon,klev) :: pfracld ! precip fraction cloudy part [-] ! OUT + REAL, DIMENSION(klon,klev) :: radocond ! condensed water used in the radiation scheme [kg/kg] ! OUT + REAL, DIMENSION(klon,klev) :: radicefrac ! ice fraction of condensed water for radiation scheme ! OUT + REAL, DIMENSION(klon,klev) :: rhcl ! clear-sky relative humidity [-] ! OUT + REAL, DIMENSION(klon) :: rain ! surface large-scale rainfall [kg/s/m2] ! OUT + REAL, DIMENSION(klon) :: snow ! surface large-scale snowfall [kg/s/m2] ! OUT + REAL, DIMENSION(klon,klev) :: qsatl ! saturation specific humidity wrt liquid [kg/kg] ! OUT + REAL, DIMENSION(klon,klev) :: qsats ! saturation specific humidity wrt ice [kg/kg] ! OUT + REAL, DIMENSION(klon,klev+1) :: prfl ! large-scale rainfall flux in the column [kg/s/m2] ! OUT + REAL, DIMENSION(klon,klev+1) :: psfl ! large-scale snowfall flux in the column [kg/s/m2] ! OUT + REAL, DIMENSION(klon,klev) :: distcltop ! distance to cloud top [m] ! OUT + REAL, DIMENSION(klon,klev) :: temp_cltop ! temperature of cloud top [K] ! OUT + REAL, DIMENSION(klon,klev) :: frac_impa ! scavenging fraction due tu impaction [-] ! OUT + REAL, DIMENSION(klon,klev) :: frac_nucl ! scavenging fraction due tu nucleation [-] ! OUT + REAL, DIMENSION(klon,klev) :: qclr ! specific total water content in clear sky region [kg/kg] ! OUT + REAL, DIMENSION(klon,klev) :: qcld ! specific total water content in cloudy region [kg/kg] ! OUT + REAL, DIMENSION(klon,klev) :: qss ! specific total water content in supersat region [kg/kg] ! OUT + REAL, DIMENSION(klon,klev) :: qvc ! specific vapor content in clouds [kg/kg] ! OUT + REAL, DIMENSION(klon,klev) :: rnebclr ! mesh fraction of clear sky [-] ! OUT + REAL, DIMENSION(klon,klev) :: rnebss ! mesh fraction of ISSR [-] ! OUT + REAL, DIMENSION(klon,klev) :: gamma_ss ! coefficient governing the ice nucleation RHi threshold [-] ! OUT + REAL, DIMENSION(klon,klev) :: Tcontr ! threshold temperature for contrail formation [K] ! OUT + REAL, DIMENSION(klon,klev) :: qcontr ! threshold humidity for contrail formation [kg/kg] ! OUT + REAL, DIMENSION(klon,klev) :: qcontr2 ! // (2nd expression more consistent with LMDZ expression of q)! OUT + REAL, DIMENSION(klon,klev) :: fcontrN ! fraction of grid favourable to non-persistent contrails ! OUT + REAL, DIMENSION(klon,klev) :: fcontrP ! fraction of grid favourable to persistent contrails ! OUT +!-------------------------------------------------------------- + + REAL, DIMENSION(klon,klev) :: d_t_eva,d_q_eva,d_ql_eva,d_qi_eva +include "YOMCST.h" ! include "clesphys.h" @@ -58,15 +162,22 @@ -real :: temp_newton(klon,klev) -integer :: k +real,dimension(klon,klev) :: temp_newton +integer :: i,k,iq +INTEGER, SAVE :: itap=0 +!$OMP THREADPRIVATE(itap) +INTEGER, SAVE :: abortphy=0 ! Reprere si on doit arreter en fin de phys +!$OMP THREADPRIVATE(abortphy) + +integer, save :: iflag_reevap=1,iflag_newton=0,iflag_vdif=1,iflag_lscp=1,iflag_cloudth_vert=3,iflag_ratqs=4 +!$OMP THREADPRIVATE(iflag_reevap,iflag_newton,iflag_vdif,iflag_lscp,iflag_cloudth_vert,iflag_ratqs) + logical, save :: first=.true. !$OMP THREADPRIVATE(first) - -real,save :: rg=9.81 -!$OMP THREADPRIVATE(rg) ! For I/Os integer :: itau0 real :: zjulian +real,dimension(klon,klev) :: du0,dv0,dqbs0 +real,dimension(klon,klev) :: cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv @@ -93,4 +204,43 @@ ! Initialize IOIPSL output file #endif +call suphel +call vdif_ini_(klon,RCPD,RD,RG,RKAPPA) +! Pourquoi ce tau_thermals en argument ??? AFAIRE +tau_thermals=0. +call getin_p('iflag_thermals',iflag_thermals) + +call getin_p('iflag_newton',iflag_newton) +call getin_p('iflag_reevap',iflag_reevap) +call getin_p('iflag_cloudth_vert',iflag_cloudth_vert) +call getin_p('iflag_ratqs',iflag_ratqs) +call getin_p('iflag_vdif',iflag_vdif) +call getin_p('iflag_lscp',iflag_lscp) +call getin_p('ratqsbas',ratqsbas) +call getin_p('ratqshaut',ratqshaut) +call getin_p('ratqsp0',ratqsp0) +call getin_p('ratqsdp',ratqsdp) +CALL thermcell_ini(iflag_thermals,0,tau_thermals,6, & + & RG,RD,RCPD,RKAPPA,RLVTT,RETV) +CALL lscp_ini(pdtphys,.false.,iflag_ratqs, RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG) + + + +allocate(tsrf(klon),q2(klon,klev+1),f0(klon),zmax0(klon)) +allocate(u_prev(klon,klev),v_prev(klon,klev),t_prev(klon,klev),q_prev(klon,klev,nqtot)) + +u_prev(:,:)=u(:,:) +v_prev(:,:)=v(:,:) +t_prev(:,:)=temp(:,:) +q_prev(:,:,:)=qx(:,:,:) + +q2=1.e-10 +tsrf=temp(:,1) +f0=0. +zmax0=0. + +iflag_replay=0 +call getin_p('iflag_replay',iflag_replay) +if ( iflag_replay >= 0 ) CALL iophys_ini(pdtphys) + endif ! of if (debut) @@ -100,4 +250,8 @@ !------------------------------------------------------------ +d_u_dyn(:,:)=(u(:,:)-u_prev(:,:))/pdtphys +d_v_dyn(:,:)=(v(:,:)-v_prev(:,:))/pdtphys +d_t_dyn(:,:)=(temp(:,:)-t_prev(:,:))/pdtphys +d_q_dyn(:,:,:)=(qx(:,:,:)-q_prev(:,:,:))/pdtphys ! set all tendencies to zero @@ -108,18 +262,280 @@ d_ps(1:klon)=0. +u_loc(1:klon,1:klev)=u(1:klon,1:klev) +v_loc(1:klon,1:klev)=v(1:klon,1:klev) +t_loc(1:klon,1:klev)=temp(1:klon,1:klev) +d_u_loc(1:klon,1:klev)=0. +d_v_loc(1:klon,1:klev)=0. +d_t_loc(1:klon,1:klev)=0. +do iq=1,nqtot + do k=1,klev + do i=1,klon + q_loc(i,k,iq)=qx(i,k,iq) + enddo + enddo +enddo + +du0(1:klon,1:klev)=0. +dv0(1:klon,1:klev)=0. +dqbs0(1:klon,1:klev)=0. + + + !------------------------------------------------------------ ! Calculs !------------------------------------------------------------ -! compute tendencies to return to the dynamics: -! "friction" on the first layer -d_u(1:klon,1)=-u(1:klon,1)/86400. -d_v(1:klon,1)=-v(1:klon,1)/86400. -! newtonian relaxation towards temp_newton() +!------------------------------------------------------------ +! Rappel en temperature et frottement dans la premiere chouche +!------------------------------------------------------------ + +if ( iflag_newton == 1 ) then + ! compute tendencies to return to the dynamics: + ! "friction" on the first layer + d_u(1:klon,1)=-u(1:klon,1)/86400. + d_v(1:klon,1)=-v(1:klon,1)/86400. + ! newtonian relaxation towards temp_newton() + do k=1,klev + temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3 + d_t(1:klon,k)=(temp_newton(1:klon,k)-temp(1:klon,k))/5.e5 + enddo +else + temp_newton(:,:)=0. +endif + + +!------------------------------------------------------------ +! Reevaporation de la pluie +!------------------------------------------------------------ + +iflag_ice_thermo=1 +if ( iflag_reevap == 1 ) then + CALL reevap (klon,klev,iflag_ice_thermo,t_loc,q_loc(:,:,1),q_loc(:,:,2),q_loc(:,:,3), & +& d_t_eva,d_q_eva,d_ql_eva,d_qi_eva) + do k=1,klev + do i=1,klon + t_loc(i,k)=t_loc(i,k)+d_t_eva(i,k) + q_loc(i,k,1)=q_loc(i,k,1)+d_q_eva(i,k) + q_loc(i,k,2)=q_loc(i,k,2)+d_ql_eva(i,k) + q_loc(i,k,3)=q_loc(i,k,3)+d_qi_eva(i,k) + q_loc(i,k,2)=0. + q_loc(i,k,3)=0. + enddo + enddo +else + d_t_eva(:,:)=0. + d_q_eva(:,:)=0. + d_ql_eva(:,:)=0. + d_qi_eva(:,:)=0. +endif + + + +!----------------------------------------------------------------------- +! Variables intermédiaires (altitudes, temperature potentielle ...) +!----------------------------------------------------------------------- + +DO k=1,klev + DO i=1,klon + zzlay(i,k)=pphi(i,k)/rg + ENDDO +ENDDO +DO i=1,klon + zzlev(i,1)=0. +ENDDO +DO k=2,klev + DO i=1,klon + z1=(pplay(i,k-1)+paprs(i,k))/(pplay(i,k-1)-paprs(i,k)) + z2=(paprs(i,k)+pplay(i,k))/(paprs(i,k)-pplay(i,k)) + zzlev(i,k)=(z1*zzlay(i,k-1)+z2*zzlay(i,k))/(z1+z2) + ENDDO +ENDDO + +! Transformation de la temperature en temperature potentielle +DO k=1,klev + DO i=1,klon + zpopsk(i,k)=(pplay(i,k)/paprs(i,1))**rkappa + masse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg + ENDDO +ENDDO +DO k=1,klev + DO i=1,klon + h_loc(i,k)=t_loc(i,k)/zpopsk(i,k) + d_h_loc(i,k)=d_t_loc(i,k)/zpopsk(i,k) + d_q_loc(i,k,1)=0. + ENDDO +ENDDO + +!----------------------------------------------------------------------- +! Diffusion verticale +!----------------------------------------------------------------------- + +if ( iflag_vdif == 1 ) then + emis(:)=1. + !tsrf=300. + z0m=0.035 + z0h=0.035 + capcal=1e2 + lwrite=.false. + print*,'lwrite ',lwrite + call vdif(klon,klev, & + & pdtphys,capcal,z0m,z0h, & + & pplay,paprs,zzlay,zzlev, & + & u_loc,v_loc,t_loc,h_loc,q_loc,tsrf,emis, & + & d_u_loc,d_v_loc,d_h_loc,d_q_loc,fluxsrf, & + & d_u_vdif,d_v_vdif,d_h_vdif,d_q_vdif,dtsrf,q2,kz_v,kz_h, & + & richardson,cdv,cdh, & + & lwrite) + do k=1,klev + do i=1,klon + d_t_vdif(i,k)=d_h_vdif(i,k)*zpopsk(i,k) + t_loc(i,k)=t_loc(i,k)+d_t_vdif(i,k)*pdtphys + u_loc(i,k)=u_loc(i,k)+d_u_vdif(i,k)*pdtphys + v_loc(i,k)=v_loc(i,k)+d_v_vdif(i,k)*pdtphys + q_loc(i,k,1)=q_loc(i,k,1)+d_q_vdif(i,k,1)*pdtphys + enddo + enddo + do i=1,klon + tsrf(i)=tsrf(i)+dtsrf(i)*pdtphys + enddo +else + d_u_vdif(:,:)=0. + d_v_vdif(:,:)=0. + d_t_vdif(:,:)=0. + d_h_vdif(:,:)=0. + d_q_vdif(:,:,1)=0. + kz_v(:,:)=0. + kz_h(:,:)=0. + richardson(:,:)=0. +endif + +!----------------------------------------------------------------------- +! Thermiques +!----------------------------------------------------------------------- + do k=1,klev - temp_newton(1:klon,k)=280.+cos(latitude(1:klon))*40.-pphi(1:klon,k)/rg*6.e-3 - d_t(1:klon,k)=(temp_newton(1:klon,k)-t(1:klon,k))/1.e5 + do i=1,klon + d_u_the(i,k)=0. + d_v_the(i,k)=0. + d_t_the(i,k)=0. + d_q_the(i,k,1)=0. + enddo enddo +if ( iflag_thermals > 0 ) then + + + zqta(:,:)=q_loc(:,:,1) + call caltherm(pdtphys & + & ,pplay,paprs,pphi & + & ,u_loc,v_loc,t_loc,q_loc,debut & + & ,f0,zmax0,d_u_the,d_v_the,d_t_the,d_q_the & + & ,frac_the,fm_the,zqta,ztv,zpspsk,ztla,zthl & + & ) + + do k=1,klev + do i=1,klon + t_loc(i,k)=t_loc(i,k)+d_t_the(i,k) + u_loc(i,k)=u_loc(i,k)+d_u_the(i,k) + v_loc(i,k)=v_loc(i,k)+d_v_the(i,k) + q_loc(i,k,1)=q_loc(i,k,1)+d_q_the(i,k,1) + enddo + enddo + +else + frac_the(:,:)=0. + fm_the(:,:)=0. + ztv(:,:)=t_loc(:,:) + zqta(:,:)=q_loc(:,:,1) + ztla(:,:)=0. + zthl(:,:)=0. + zpspsk(:,:)=(pplay(:,:)/100000.)**RKAPPA + +endif + +!----------------------------------------------------------------------- +! Condensation grande échelle +!----------------------------------------------------------------------- + +iflag_cld_th=5 +ok_ice_sursat=.false. +ptconv(:,:)=.false. +distcltop=0. +temp_cltop=0. +beta(:,:)=1. +rneb_seri(:,:)=0. +do k=1,klev + ratqs(:,k)=ratqsbas+0.5*(ratqshaut-ratqsbas) & + *( tanh( (ratqsp0-pplay(:,k))/ratqsdp) + 1.) +enddo + + +if ( iflag_lscp == 1 ) then + + call lscp(klon,klev,pdtphys,missing_val, & + paprs,pplay,t_loc,q_loc,ptconv,ratqs, & + d_t_lscp, d_q_lscp, d_ql_lscp, d_qi_lscp, rneb, rneblsvol, rneb_seri, & + pfraclr,pfracld, & + radocond, radicefrac, rain, snow, & + frac_impa, frac_nucl, beta, & + prfl, psfl, rhcl, zqta, frac_the, & + ztv, zpspsk, ztla, zthl, iflag_cld_th, & + iflag_ice_thermo, ok_ice_sursat, qsatl, qsats, & + distcltop,temp_cltop, & + qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, & + Tcontr, qcontr, qcontr2, fcontrN, fcontrP, & + cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv) + + + do k=1,klev + do i=1,klon + t_loc(i,k)=t_loc(i,k)+d_t_lscp(i,k) + q_loc(i,k,1)=q_loc(i,k,1)+d_q_lscp(i,k) + q_loc(i,k,2)=q_loc(i,k,2)+d_ql_lscp(i,k) + q_loc(i,k,3)=q_loc(i,k,3)+d_qi_lscp(i,k) + enddo + enddo + +else + d_t_lscp(:,:)=0. + d_q_lscp(:,:)=0. + d_ql_lscp(:,:)=0. + d_qi_lscp(:,:)=0. + rneb(:,:)=0. + rneblsvol(:,:)=0. + pfraclr(:,:)=0. + pfracld(:,:)=0. + radocond(:,:)=0. + rain(:)=0. + snow(:)=0. + radicefrac(:,:)=0. + rhcl (:,:)=0. + qsatl (:,:)=0. + qsats (:,:)=0. + prfl (:,:)=0. + psfl (:,:)=0. + distcltop (:,:)=0. + temp_cltop(:,:)=0. + frac_impa (:,:)=0. + frac_nucl (:,:)=0. + qclr (:,:)=0. + qcld (:,:)=0. + qss (:,:)=0. + qvc (:,:)=0. + rnebclr (:,:)=0. + rnebss (:,:)=0. + gamma_ss (:,:)=0. + Tcontr (:,:)=0. + qcontr (:,:)=0. + qcontr2 (:,:)=0. + fcontrN (:,:)=0. + fcontrP (:,:)=0. +endif + + +d_u(:,:)=(u_loc(:,:)-u(:,:))/pdtphys +d_v(:,:)=(v_loc(:,:)-v(:,:))/pdtphys +d_t(:,:)=(t_loc(:,:)-temp(:,:))/pdtphys +d_qx(:,:,:)=(q_loc(:,:,:)-qx(:,:,:))/pdtphys !------------------------------------------------------------ @@ -128,5 +544,15 @@ -call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,t,qx,0.*t,0.*t,0.*t,0.*t,0.*t,0.*t) +tsrf_(:)=tsrf(:) +if ( iflag_replay == -1 ) then + call output_physiqex(debut,zjulian,pdtphys,presnivs,paprs,u,v,temp,qx,0.*u,0.*u,0.*u,0.*u,q2,0.*u) +else if (iflag_replay == 0 ) then + ! En mode replay, on sort aussi les variables de base +! Les lignes qui suivent ont été générées automatiquement avec : +! ( for i in `grep -i 'real.*::' physiqex_mod.F90 | sed -e '/^!/d' | grep '(klon,klev' | cut -d: -f3 | cut -d! -f1 | sed -e 's/,/ /g' -e '/rot/d'` ; do echo ' call iophys_ecrit("'$i'",klev,"","",'$i')' ; done ) > physiqex_out.h +! ( for i in `grep -i 'real.*::' physiqex_mod.F90 | sed -e '/^!/d' | grep '(klon)' | cut -d: -f3 | cut -d! -f1 | sed -e 's/,/ /g' -e '/rot/d'` ; do echo ' call iophys_ecrit("'$i'",1,"","",'$i')' ; done ) >> physiqex_out.h +include "physiqex_out.h" + +endif @@ -136,4 +562,5 @@ endif +print*,'Fin physiqex' end subroutine physiqex Index: LMDZ6/trunk/libf/phylmd/vdif_kcay.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/vdif_kcay.F90 (revision 4653) +++ LMDZ6/trunk/libf/phylmd/vdif_kcay.F90 (revision 4654) @@ -2,7 +2,7 @@ ! $Header$ -SUBROUTINE vdif_kcay(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, & +SUBROUTINE vdif_kcay(klon,klev, dt, g, rconst, plev, temp, zlev, zlay, u, v, & teta, cd, q2, q2diag, km, kn, ustar, l_mix) - USE dimphy + IMPLICIT NONE @@ -28,20 +28,26 @@ ! ....................................................................... - REAL dt, g, rconst - REAL plev(klon, klev+1), temp(klon, klev) - REAL ustar(klon), snstable - REAL zlev(klon, klev+1) - REAL zlay(klon, klev) - REAL u(klon, klev) - REAL v(klon, klev) - REAL teta(klon, klev) - REAL cd(klon) - REAL q2(klon, klev+1), q2s(klon, klev+1) - REAL q2diag(klon, klev+1) - REAL km(klon, klev+1) - REAL kn(klon, klev+1) - REAL sq(klon), sqz(klon), zz(klon, klev+1), zq, long0(klon) - - INTEGER l_mix, iii + INTEGER, INTENT(IN) :: klon,klev + REAL,INTENT(IN) :: dt, g, rconst + REAL,DIMENSION(klon,klev+1),INTENT(IN) :: plev + REAL,DIMENSION(klon,klev),INTENT(IN) :: temp + REAL,DIMENSION(klon),INTENT(IN) :: ustar(klon) + REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: zlev + REAL,DIMENSION(klon,klev),INTENT(IN) :: zlay + REAL,DIMENSION(klon,klev),INTENT(IN) :: u + REAL,DIMENSION(klon,klev),INTENT(IN) :: v + REAL,DIMENSION(klon,klev),INTENT(IN) :: teta + REAL,DIMENSION(klon),INTENT(IN) :: cd(klon) + REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: q2 + REAL,DIMENSION(klon,klev+1),INTENT(INOUT) :: q2diag + REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: km + REAL,DIMENSION(klon,klev+1),INTENT(OUT) :: kn + INTEGER, INTENT(OUT) :: l_mix + + REAL,DIMENSION(klon) :: sq, sqz, long0 + REAL,DIMENSION(klon,klev+1) :: q2s,zz + REAL :: snstable,zq + + INTEGER :: iii ! ....................................................................... @@ -56,8 +62,7 @@ ! ....................................................................... - INTEGER nlay, nlev, ngrid - REAL unsdz(klon, klev) - REAL unsdzdec(klon, klev+1) - REAL q(klon, klev+1) + INTEGER :: nlay, nlev, ngrid + REAL, DIMENSION(klon,klev) :: unsdz + REAL, DIMENSION(klon, klev+1) :: unsdzdec,q ! ....................................................................... @@ -70,8 +75,8 @@ ! ....................................................................... - REAL kmpre(klon, klev+1) - REAL qcstat - REAL q2cstat - REAL sss, sssq + REAL, DIMENSION(klon, klev+1) :: kmpre + REAL :: qcstat + REAL :: q2cstat + REAL :: sss, sssq ! ....................................................................... @@ -79,5 +84,5 @@ ! ....................................................................... - REAL long(klon, klev+1) + REAL, DIMENSION(klon, klev+1) :: long ! ....................................................................... @@ -97,13 +102,10 @@ ! ....................................................................... - REAL kmq3 - REAL kmcstat - REAL knq3 - REAL mcstat - REAL m2cstat - REAL m(klon, klev+1) - REAL mpre(klon, klev+1) - REAL m2(klon, klev+1) - REAL n2(klon, klev+1) + REAL :: kmq3 + REAL :: kmcstat + REAL :: knq3 + REAL :: mcstat + REAL :: m2cstat + REAL, DIMENSION(klon, klev+1) :: m,mpre,m2,n2 ! ....................................................................... @@ -121,15 +123,10 @@ ! ....................................................................... - REAL gn - REAL gnmin - REAL gnmax - LOGICAL gninf - LOGICAL gnsup - REAL gm - ! REAL ri(klon,klev+1) - REAL sn(klon, klev+1) - REAL snq2(klon, klev+1) - REAL sm(klon, klev+1) - REAL smq2(klon, klev+1) + REAL :: gn,gm + REAL :: gnmin + REAL :: gnmax + LOGICAL :: gninf + LOGICAL :: gnsup + REAL, DIMENSION(klon, klev+1) :: sn, snq2, sm, smq2 ! ....................................................................... @@ -142,9 +139,9 @@ ! ....................................................................... - REAL kappa - REAL long00 - REAL a1, a2, b1, b2, c1 - REAL cn1, cn2 - REAL cm1, cm2, cm3, cm4 + REAL :: kappa + REAL :: long00 + REAL :: a1, a2, b1, b2, c1 + REAL :: cn1, cn2 + REAL :: cm1, cm2, cm3, cm4 ! ....................................................................... @@ -155,8 +152,8 @@ ! ....................................................................... - REAL termq - REAL termq3 - REAL termqm2 - REAL termq3m2 + REAL :: termq + REAL :: termq3 + REAL :: termqm2 + REAL :: termq3m2 ! ....................................................................... @@ -165,15 +162,15 @@ ! ....................................................................... - REAL q2min - REAL q2max + REAL :: q2min + REAL :: q2max ! ....................................................................... ! knmin : borne inferieure de kn ! kmmin : borne inferieure de km ! ....................................................................... - REAL knmin - REAL kmmin - ! ....................................................................... - INTEGER ilay, ilev, igrid - REAL tmp1, tmp2 + REAL :: knmin + REAL :: kmmin + ! ....................................................................... + INTEGER :: ilay, ilev, igrid + REAL :: tmp1, tmp2 ! ....................................................................... PARAMETER (kappa=0.4E+0) @@ -202,5 +199,5 @@ PARAMETER (cm4=-9.E+0*a1*a2) - LOGICAL first + LOGICAL :: first SAVE first DATA first/.TRUE./ @@ -211,13 +208,16 @@ ! Initialisation de q2 +ngrid=klon nlay = klev nlev = klev + 1 - CALL yamada(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, cd, & - q2diag, km, kn, ustar, l_mix) - IF (first .AND. 1==1) THEN - first = .FALSE. - q2 = q2diag - END IF +! Initialisation avec un schema d'equilibre +! CALL yamada(ngrid, dt, g, rconst, plev, temp, zlev, zlay, u, v, teta, cd, & +! q2diag, km, kn, ustar, l_mix) +! IF (first .AND. 1==1) THEN +! first = .FALSE. +! q2 = q2diag +! END IF +q2diag=0. DO ilev = 2, nlay