Index: LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90 (revision 4558) +++ LMDZ6/trunk/libf/phylmd/lscp_ini_mod.F90 (revision 4559) @@ -12,6 +12,6 @@ !$OMP THREADPRIVATE(seuil_neb) - INTEGER, SAVE :: ninter=5 ! number of iterations to calculate autoconversion to precipitation - !$OMP THREADPRIVATE(ninter) + INTEGER, SAVE :: niter_lscp=5 ! number of iterations to calculate autoconversion to precipitation + !$OMP THREADPRIVATE(niter_lscp) INTEGER,SAVE :: iflag_evap_prec=1 ! precipitation evaporation flag. 0: nothing, 1: "old way", @@ -41,7 +41,4 @@ !$OMP THREADPRIVATE(ok_radocond_snow) - LOGICAL, SAVE :: ok_debug_autoconversion=.true. ! removes a bug in the autoconversion process - !$OMP THREADPRIVATE(ok_debug_autoconversion) - REAL, SAVE :: t_glace_min=258.0 ! lower-bound temperature parameter for cloud phase determination !$OMP THREADPRIVATE(t_glace_min) @@ -77,4 +74,7 @@ !$OMP THREADPRIVATE(iflag_pdf) + INTEGER, SAVE :: iflag_autoconversion=0 ! autoconversion option + !$OMP THREADPRIVATE(iflag_autoconversion) + LOGICAL, SAVE :: reevap_ice=.false. ! no liquid precip for T< threshold !$OMP THREADPRIVATE(reevap_ice) @@ -92,4 +92,10 @@ !$OMP THREADPRIVATE(cld_tau_con) + REAL, SAVE :: cld_expo_lsc=2. ! liquid autoconversion threshold exponent, stratiform rain + !$OMP THREADPRIVATE(cld_expo_lsc) + + REAL, SAVE :: cld_expo_con=2. ! liquid autoconversion threshold exponent, convective rain + !$OMP THREADPRIVATE(cld_expo_con) + REAL, SAVE :: ffallv_lsc=1. ! tuning coefficient crystal fall velocity, stratiform !$OMP THREADPRIVATE(ffallv_lsc) @@ -104,5 +110,9 @@ !$OMP THREADPRIVATE(coef_eva_i) - + REAL cice_velo ! factor in the ice fall velocity formulation + PARAMETER (cice_velo=1.645) + + REAL dice_velo ! exponent in the ice fall velocity formulation + PARAMETER (dice_velo=0.16) @@ -122,4 +132,6 @@ REAL, INTENT(IN) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in REAL, INTENT(IN) :: RVTMP2_in, RTT_in, RD_in, RG_in + character (len=20) :: modname='lscp_ini_mod' + character (len=80) :: abort_message @@ -136,5 +148,5 @@ - CALL getin_p('ninter',ninter) + CALL getin_p('niter_lscp',niter_lscp) CALL getin_p('iflag_evap_prec',iflag_evap_prec) CALL getin_p('seuil_neb',seuil_neb) @@ -142,5 +154,4 @@ CALL getin_p('iflag_mpc_bl',iflag_mpc_bl) CALL getin_p('ok_radocond_snow',ok_radocond_snow) - CALL getin_p('ok_debug_autoconversion',ok_debug_autoconversion) CALL getin_p('t_glace_max',t_glace_max) CALL getin_p('t_glace_min',t_glace_min) @@ -159,4 +170,6 @@ CALL getin_p('cld_tau_lsc',cld_tau_lsc) CALL getin_p('cld_tau_con',cld_tau_con) + CALL getin_p('cld_expo_lsc',cld_expo_lsc) + CALL getin_p('cld_expo_con',cld_expo_con) CALL getin_p('ffallv_lsc',ffallv_lsc) CALL getin_p('ffallv_lsc',ffallv_con) @@ -164,9 +177,10 @@ coef_eva_i=coef_eva CALL getin_p('coef_eva_i',coef_eva_i) - - - - - WRITE(lunout,*) 'lscp, ninter:', ninter + CALL getin_p('iflag_autoconversion',iflag_autoconversion) + + + + + WRITE(lunout,*) 'lscp, niter_lscp:', niter_lscp WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb @@ -174,5 +188,4 @@ WRITE(lunout,*) 'lscp, iflag_mpc_bl:', iflag_mpc_bl WRITE(lunout,*) 'lscp, ok_radocond_snow:', ok_radocond_snow - WRITE(lunout,*) 'lscp, ok_debug_autoconversion:', ok_debug_autoconversion WRITE(lunout,*) 'lscp, t_glace_max:', t_glace_max WRITE(lunout,*) 'lscp, t_glace_min:', t_glace_min @@ -191,8 +204,11 @@ WRITE(lunout,*) 'lscp, cld_tau_lsc', cld_tau_lsc WRITE(lunout,*) 'lscp, cld_tau_con', cld_tau_con + WRITE(lunout,*) 'lscp, cld_expo_lsc', cld_expo_lsc + WRITE(lunout,*) 'lscp, cld_expo_con', cld_expo_con WRITE(lunout,*) 'lscp, ffallv_lsc', ffallv_lsc WRITE(lunout,*) 'lscp, ffallv_con', ffallv_con WRITE(lunout,*) 'lscp, coef_eva', coef_eva WRITE(lunout,*) 'lscp, coef_eva_i', coef_eva_i + WRITE(lunout,*) 'lscp, iflag_autoconversion', iflag_autoconversion @@ -201,7 +217,17 @@ ! check for precipitation sub-time steps - IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN + IF (ABS(dtime/REAL(niter_lscp)-360.0).GT.0.001) THEN WRITE(lunout,*) 'lscp: it is not expected, see Z.X.Li', dtime WRITE(lunout,*) 'I would prefer a 6 min sub-timestep' + ENDIF + + ! check consistency between numerical resolution of autoconversion + ! and other options + + IF (iflag_autoconversion .EQ. 2) THEN + IF ((iflag_vice .NE. 0) .OR. (niter_lscp .GT. 1)) THEN + abort_message = 'in lscp, iflag_autoconversion=2 requires iflag_vice=0 and niter_lscp=1' + CALL abort_physic (modname,abort_message,1) + ENDIF ENDIF Index: LMDZ6/trunk/libf/phylmd/lscp_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lscp_mod.F90 (revision 4558) +++ LMDZ6/trunk/libf/phylmd/lscp_mod.F90 (revision 4559) @@ -93,9 +93,10 @@ USE ice_sursat_mod, ONLY : ice_sursat -USE lscp_ini_mod, ONLY : seuil_neb, ninter, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min -USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, ok_debug_autoconversion +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 USE lscp_ini_mod, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min USE lscp_ini_mod, ONLY : coef_eva, coef_eva_i,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con -USE lscp_ini_mod, ONLY : iflag_bergeron, iflag_fisrtilp_qsat +USE lscp_ini_mod, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo +USE lscp_ini_mod, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc USE lscp_ini_mod, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG @@ -190,5 +191,5 @@ REAL qsl(klon), qsi(klon) - REAL zct, zcl + REAL zct, zcl,zexpo REAL ctot(klon,klev) REAL ctot_vol(klon,klev) @@ -232,7 +233,8 @@ 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 + REAL velo(klon,klev), vr, ffallv REAL qlmpc, qimpc, rnebmpc REAL radocondi(klon,klev), radocondl(klon,klev) + REAL effective_zneb INTEGER i, k, n, kk @@ -922,4 +924,9 @@ ENDIF + + ! Autoconversion + ! ------------------------------------------------------------------------------- + + ! Initialisation of zoliq and radocond variables @@ -933,14 +940,12 @@ iwc(i) = 0. zneb(i) = MAX(rneb(i,k),seuil_neb) - radocond(i,k) = zoliq(i)/REAL(ninter+1) + radocond(i,k) = zoliq(i)/REAL(niter_lscp+1) radicefrac(i,k) = zfice(i) - radocondi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1) - radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1) + radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1) + radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1) ENDDO - ! Autoconversion - ! ------------------------------------------------------------------------------- - - DO n = 1, ninter + + DO n = 1, niter_lscp DO i=1,klon @@ -968,18 +973,16 @@ ELSE - IF (ok_debug_autoconversion) THEN ! if condition to be removed after test phase - - ! water quantity to remove: zchau (Sundqvist, 1978) - ! same thing for the ice: zfroi (Zender & Kiehl, 1997) IF (ptconv(i,k)) THEN ! if convective point zcl=cld_lc_con zct=1./cld_tau_con - zfroi = dtime/REAL(ninter)/zdz(i)*zoliqi(i)*velo(i,k) + zexpo=cld_expo_con + ffallv=ffallv_con ELSE zcl=cld_lc_lsc zct=1./cld_tau_lsc - zfroi = dtime/REAL(ninter)/zdz(i)*zoliqi(i) & ! dqice/dt=1/rho*d(rho*wice*qice)/dz - *velo(i,k) + zexpo=cld_expo_lsc + ffallv=ffallv_lsc ENDIF + ! if vertical heterogeneity is taken into account, we use @@ -988,26 +991,8 @@ ! reduces the in-cloud water). + ! Liquid water quantity to remove: zchau (Sundqvist, 1978) + ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2) + !......................................................... IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN - zchau = zct *dtime/REAL(ninter) * zoliql(i) & - *(1.0-EXP(-(zoliql(i)/ctot_vol(i,k)/zcl)**2)) - ELSE - zchau = zct *dtime/REAL(ninter) * zoliql(i) & - *(1.0-EXP(-(zoliql(i)/zneb(i)/zcl)**2)) ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2) - ENDIF - - ELSE ! with old bug in autoconversion - - ! water quantity to remove: zchau (Sundqvist, 1978) - ! same thing for the ice: zfroi (Zender & Kiehl, 1997) - IF (ptconv(i,k)) THEN ! if convective point - zcl=cld_lc_con - zct=1./cld_tau_con - zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i,k)*zfice(i) - ELSE - zcl=cld_lc_lsc - zct=1./cld_tau_lsc - zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) & ! dqice/dt=1/rho*d(rho*wice*qice)/dz - *velo(i,k) * zfice(i) - ENDIF ! if vertical heterogeneity is taken into account, we use @@ -1015,15 +1000,39 @@ ! surface fraction (which is larger and artificially ! reduces the in-cloud water). - - - IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN - zchau = zct *dtime/REAL(ninter) * zoliq(i) & - *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl)**2)) *(1.-zfice(i)) + effective_zneb=ctot_vol(i,k) ELSE - zchau = zct *dtime/REAL(ninter) * zoliq(i) & - *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl)**2)) *(1.-zfice(i)) ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2) + effective_zneb=zneb(i) + ENDIF + + + IF (iflag_autoconversion .EQ. 2) THEN + ! two-steps resolution with niter_lscp=1 sufficient + ! we first treat the second term (with exponential) in an explicit way + ! and then treat the first term (-q/tau) in an exact way + zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct & + *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo)))) + ELSE + ! old explicit resolution with subtimesteps + zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) & + *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo)) ENDIF - ENDIF ! ok_debug_autoconversion + + ! Ice water quantity to remove (Zender & Kiehl, 1997) + ! dqice/dt=1/rho*d(rho*wice*qice)/dz + !.................................... + IF (iflag_autoconversion .EQ. 2) THEN + ! exact resolution, niter_lscp=1 is sufficient but works only + ! with iflag_vice=0 + IF (zoliqi(i) .GT. 0.) THEN + zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) & + +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo)) + ELSE + zfroi=0. + ENDIF + ELSE + ! old explicit resolution with subtimesteps + zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*velo(i,k) + ENDIF zrain = MIN(MAX(zchau,0.0),zoliql(i)) @@ -1033,13 +1042,21 @@ ENDIF - zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain , 0.0) - zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow , 0.0) - zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) + IF (iflag_autoconversion .GE. 1) THEN + ! debugged version with phase conservation through the autoconversion process + zoliql(i) = MAX(zoliql(i)-1.*zrain , 0.0) + zoliqi(i) = MAX(zoliqi(i)-1.*zsnow , 0.0) + zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) + ELSE + ! bugged version with phase resetting + zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain , 0.0) + zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow , 0.0) + zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) + ENDIF ! c_iso: call isotope_conversion (for readibility) or duplicate - radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(ninter+1) - radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(ninter+1) - radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(ninter+1) + radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1) + radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1) + radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1) ENDIF ! rneb >0 @@ -1137,4 +1154,5 @@ ENDDO + ENDIF Index: LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90 (revision 4558) +++ LMDZ6/trunk/libf/phylmd/lscp_tools_mod.F90 (revision 4559) @@ -17,4 +17,5 @@ use lscp_ini_mod, only: iflag_vice, ffallv_con, ffallv_lsc + use lscp_ini_mod, only: cice_velo, dice_velo IMPLICIT NONE @@ -31,5 +32,5 @@ INTEGER i - REAL logvm,iwcg,tempc,phpa,cvel,dvel,fallv_tun + REAL logvm,iwcg,tempc,phpa,fallv_tun REAL m2ice, m2snow, vmice, vmsnow REAL aice, bice, asnow, bsnow @@ -60,6 +61,4 @@ velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s - dvel=0.2 - cvel=fallv_tun*65.0*(rho(i)**0.2)*(150./phpa)**0.15 ELSE IF (iflag_vice .EQ. 2) THEN @@ -92,12 +91,8 @@ vmsnow=vmsnow*((1000./phpa)**0.35) velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s - dvel=0.2 - cvel=velo(i)/((iwcg/1000.*rho(i))**dvel) ELSE ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990 - velo(i) = fallv_tun*3.29/2.0*((iwcg/1000.)**0.16) - dvel=0.16 - cvel=fallv_tun*3.29/2.0*(rho(i)**0.16) + velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo) ENDIF ENDDO